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ABSTRACT 

Two methods are developed for solving the steady-state spherically symmetric radia- 
tive transfer equation for resonance line radiation emitted by a point source in the 
Intergalactic Medium, in the context of the Wouthuysen-Field mechanism for cou- 

■ pling the hyperfine structure spin temperature of hydrogen to the gas temperature. 

One method is based on solving the ray and moment equations using finite differences. 
The second uses a Monte Carlo approach incorporating methods that greatly improve 

■ the accuracy compared with previous approaches in this context. Several applications 

are presented serving as test problems for both a static medium and an expanding 

q medium, including inhomogeneities in the density and velocity fields. Solutions are 

'_- ■ obtained in the coherent scattering limit and for Doppler RII redistribution with and 

| without recoils. We find generally that the radiation intensity is linear in the cosine 

of the azimuthal angle with respect to radius to high accuracy over a broad frequency 
region across the line centre for both linear and perturbed velocity fields, yielding 
the Eddington factors /„ ~ 1/3 and g v ~ 3/5. The radiation field produced by a 
point source divides into three spatial regimes for a uniformly expanding homoge- 
' neous medium. The regimes are governed by the fraction of the distance r from the 

(/*■) . source in terms of the distance r* required for a photon to redshift from line centre 

f*** ■ to the frequency needed to escape from the expanding gas. For a standard cosmology, 

^sO | before the Universe was reionized r* takes on the universal value independent of red- 

. shift of 1.1 Mpc, depending only on the ratio of the baryon to dark matter density. At 

r/r* < 1, the radiation field is accurately described in the diffusion approximation, 
fN| . with the scattering rate declining with the distance from the source as r~ 7 / 3 , except 

at r/r* <C 1 where frequency redistribution nearly doubles the mean intensity around 
line centre. At r/r* > 1, the diffusion approximation breaks down and the decline of 
the mean intensity near line centre and the scattering rate approach the geometric 
' dilution scaling 1/r 2 . The mean intensity and scattering rate are found to be very sen- 

sitive to the gradient of the velocity field, growing exponentially with the amplitude 
of the perturbation as the limit of a vanishing velocity gradient is approached near 
the source. We expect the 21cm signal from the Epoch of Reionization to thus be a 
sensitive probe of both the density and the peculiar velocity fields. 

The solutions for the mean intensity are made available in machine-readable for- 
mat. 

Key words: atomic processess - cosmology: theory - line: formation - radiative 
transfer - radio lines: general - scattering 



1 INTRODUCTION 

The nature of formation of the first radiating objects in 
the Universe is one of the paramount unsolved problems in 
cosmological structure formation. Searches for the earliest 
galaxies have broken the spectroscopically confirmed red- 
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shift barrier of z = 7 |Vanzella et al.ll201ll : lOno et alj|2012l : 
ISchenker et"afl l2012i h with plausible candidates identified 
photometrically up to z < 9 ( McLurc et al. 201 lj), and pos- 
sibly as high as z ~ 10 (|Bouwens et aljfeoill ). These systems 
may well reside within the Epoch of Reionization (EoR) of 
hydrogen, as Cosmic Microwave Background (CMB) mea- 
surements suggest the e poch, if a sudden ev ent, occurred 
at z r = 10.4 ± 1.2 (lcr) (jKomatsu et al.ll201l[ ). If so, 21cm 
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emission from the diffuse fntergalactic Medium (IGM) would 
become visible near this same epoch, as the same sources 
would provide sufficient UV continuum to e xcite the line 
through the Wouythuysen-Field effec t (WFE) |Wouthuvsenl 
1 19521 ; iFieldl 1 19581 ; Wadau et al l Il997h . 

The prospect of discovering the EoR through the associ- 
ated 21cm signature from the diffuse IGM has inspired the 
development of a new generation of radio telescopes, such 
as the LOw Frequency Array (LOFARj], upgrades to the 
Giant Metrewave Radio Telescope (GMRTjfi the Murchi- 
son Widefield Array (MWA)B, the Primeval Structure Tele- 
scope/21 Centimeter Array (PaST/21CMA) 0, the Precision 
Array to Probe EoR (PAPER jf], and a possible Square Kilo- 
metre Array (SKA fl A recent review of this rapidly growing 
area is provided bv lPritchard fc Loebl (|201ll ). 

The interpretation of the signal will require modelling 
the radiative transfer of the Lyman resonance line photons, 
primarily Lya, that drive the WFE. Most estimates have 
presumed a homogeneous expanding medium. Early mod- 
elling neglected the effects of atomic recoil and of spatial 
diffusion about the emitting sources, assuming the sources 
were homogeneo usly and isotropically distribute d through- 
out the Universe (|Fieldlll959al ; lMadau et al.lll997T ). Allowing 
for recoil somewhat suppresses the Lya photon scattering 
rate by an amount dependin g on the local temperature and 
expansion rate of the IGM Ichen fc Miralda-Escudl |2004| ; 
iFurlanetto fc Pritch ard 2006). Using Monte Carlo solutions 
to the radiative transfer equation to include spatial diffu- 
sion, it was found th at near an emitting source the scattering 
rate varies as r~ 7//3 (|Chuzhov fc Zheng 2007; Scmclin et al.l 
I2007T ). more rapidly than the geometric dilution factor r 
predicted without spatial diffusion. It will be shown below 
that the steeper dependence arises generally over all dis- 
tances for which the radiative transfer of the Lya photons 
may be treated in the diffusion approximation. 

In reality the IGM is clumpy, with structures break- 
ing away from the cosmological expansion. Simple analytic 
estimates suggest that the scattering rate will be substan- 
tially modified not only by density and temperature fluc- 
tuations, but by gradients in the velocity field of the gas 
around individual sources l|Higgins fc Meiksin] |2009| ). More 
sophisticated methods are required for accurate solutions 
of the spatial and frequency dependent radiative trans- 
fer equat ion. Monte Carlo codes were developed for this 
purpose ilZheng fc Miralda-Escudel |2002|; iTasitsiomil 120061 ; 



IChuzhovfe Zheng||2007l ; ISemelin et al.ll2007h . The method 
has recently been applied to estima t e the expected cosmolog - 
ical 21cm signal l|Baek et al.ll2009l ; IVonlanthen et al.ll20lj ). 
A difficulty with the Monte Carlo technique is the limited 
resolution imposed by the restricted number of photon pack- 
ets that may be practically followed. An alternative grid- 
based method for solving the spherically-symmetric radia- 
tive transfer equat ion in the Eddingt on approximation has 
been developed by I Roy et all j|2009bh for linear flow fields. 
Scattering in the deuterium Lya resonance and the ad- 



1 www.lofar.org 

2 gmrt.ncra.tifr.res.in 



www 



haystack.mit.edu/ast / arrays / mwa 



4 web.phys.cmu.edu/~past 

5 astro.berkeley.edu/~dbacker / eor 
6 



www.skatclescopc.org 



dition of Lya photons produced in radiative cascades fol- 
lowing the scattering of higher order Lyman resonance line 
photons will modify the Lya mean intensity and scattering 
rate near a source. As these are smaller, sec ondary effects 
l|Chuzhov fc Zhendl2007l ; ISemelin et aLll2007h . they are not 
included in this paper. 

A further complication is the time required to establish 
a steady-state radiation field. The Lya photon scattering 
rate scales like t a ~ 2>.2n^T x ^ 2 s for a gas with neutral hy- 
drogen density nm and temperature T. Time-dependent ra- 
diative transfer computations suggest it may take 10 4 — 10 12 
scatterings to establish a steady-state radiation field, de- 
pending on the internal structure of the scattering system, 
including internal velo city gradients (jHiggins fc Meiksin! 
120091 ; iRov et alj [200981 ). Timescales of 10 6 - 10 9 yr, com- 
parable to or longer than the evolutionary timescale of star- 
bursts and quasars, may be required for structures that have 
broken away from the cosmic expansion, particularly in the 
presence of substantial x-ray heating. 

In this paper, we present two algorithms for solving the 
radiative transfer equation of an inhomogeneous medium, 
one based on finite-differences on a grid and the second 
Monte Carlo based. The grid-based method combines so- 
lutions to the ray equation and t he moment equat ion based 
on the method of iMihalas et all <| 19751 . 1 19761 . 1 1977ft . The in- 
clusion of solutions to the ray equation allows the sequence 
of moment equations to be closed without imposing the Ed- 
dington approximation. Since the method was developed for 
stellar atmospheres, modifications to the approach are de- 
scribed necessary to adapt the method to the problem of a 
source in the IGM. The method has a few restrictions: 1. it 
is implemented assuming spherical symmetry, 2. it assumes 
a steady state, and 3. it requires the velocity field around 
the source to be monotonically increasing or decreasing. In 
practice the latter is not a severe restriction because the pe- 
culiar velocity field only modulates the Hubble flow except 
within the turnaround radius very near the source. The as- 
sumption of spherical symmetry may not be very restrictive 
either, since the photons do not diffuse very far compared 
with the coherence length of cosmological structures. A so- 
lution along a ray in a 3D computation may therefore not 
differ much from that assuming the medium is isotropic with 
the radial properties of the ray. We are not able to test this, 
however, without a fully 3D solution to the radiative trans- 
fer equation, which is beyond the scope of this paper. As 
noted above, situations may arise in which the steady-state 
approximation will break down. 

The Monte Carlo method is similar to existing algo- 
rithms for resonance line radiation, but with two improve- 
ments. It incorporates RII frequency redistribution by in- 
teropolating on the RII redistribution function rather than 
directly computing collisions with atoms. This improves the 
speed of the computations by a factor of a few, but at the 
cost of requiring a frequency grid tailored to the particular 
problem. The second improvement is to compile the specific 
mean intensity based on the path lengths traversed by the 
photon packets rather than on frequency and position bin 
crossings. This improvement is general and substantially re- 
duces the noise in the specific mean intensity for a fixed 
number of photon packets. The method also serves as an in- 
dependent check on the solutions obtained through the ray 
and moments method. 
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The paper is organised as follows. The basic framework 
used to solve the radiative transfer equation in spherical 
symmetry using the ray and moments method on a grid is 
presented in the next section. In Section[3l we summarise the 
Monte Carlo method developed. We present the results of 
validation tests of both methods against analytic solutions 
for a homogeneous medium in Section [4] Both methods are 
applied to scattering problems in an inhomogeneous medium 
in Section [S] A summary and conclusions follow. Details of 
the source term used for the ray and moment solution are 
provided in Appendix [X] The Monte Carlo method is de- 
scribed in Appendix [B] The solution to the radiative trans- 
fer problem for a point continuum source in a uniformly 
expanding homogeneous medium in the diffusion approxi- 
mation is derived in Appendix [C] Tables of solutions to a 
test suite of problems using the ray and moments method are 
provided online in machine-readable format, as summarised 
in Appendix [Dl 



quantity. In solving the moment equations we assume an- 
gular information in the form of variable Eddington fac- 
tors in order to close the system of equations, however we 
are able to treat the dependence of the source function on 
the radiation field exactly. The iteration procedure is im- 
plemented by first seeking a solution to the ray equations 
utilising an estimate of the source function, before evalu- 
ating the variable Eddington factors from the ray solution 
and using them to obtain a solution to the moment equa- 
tions. It may then be necessary to improve on the estimate 
of the Eddington factors by repeating the ray solution steps 
with an improved source function estimate derived from the 
moment solution, a process that may be iterated until con- 
vergence. At each step the solution of the ray/moment equa- 
tions is obtained from finite-difference forms which are writ- 
ten _a£jna±rixequations. The matrix equations are provided 
by iHigginsl l|2012h . We next outline the method in greater 
detail. 



2 NUMERICAL RADIATIVE TRANSFER IN 
SPHERICAL SYMMETRY I: RAY AND 
MOMENT EQUATIONS 

2.1 Methodology 

We developed a grid-based method for obtaining precise nu- 
merical solutions to the radiative trans fer equation in spher- 
ical symme try based on the method of iMihalas et al.l (| 19751 . 
Il976l . Il977h . This method was developed for stellar atmo- 
spheres and consequently we have had to formulate the 
boundary conditions for a non-blackbody source. Here we 
provide a summary of our implementation. 

The radiative transfer equation for the specific intensity 
I v in the comoving frame for a spherically symmetric scat- 
tering medium with specific volume emissivity r) v , inverse 
attenuation len gth Xu and radial velocity V(r) is given by 
|Mihalasill97Sl ) 

»-dF + ~ a{r) [ X - M + " P{r)] ^7 (i) 



where a(r) and /3(r) are given by 
a(r) = (n>/c)^, 



p{r) 



&\nV{r) 



(2) 
(3) 



dlnr 

and fi is the cosine of the inclination relative to the radial di- 
rection of the radiation field described by I v (r,n), and fo is 
the resonance line frequency. We limit our study to expand- 
ing velocity fields V(r) satisfying V(r) > and V'(r) > 
for all r. (The method also may be applied to a monoton- 
ically decreasing radial velocity field. Non-monotonic fields 
would require combining solutions from piecewise monotonic 
fields.) The method requires solving a system of angular mo- 
ment equations derived from equation Q, together with a 
solution of equation {TJ itself along a specific set of rays. 
The approach is to solve the ray and moment equations si- 
multaneously, iterating between the two until adequate con- 
vergence is achieved. In solving the ray equations, we treat 
the angular dependence of the radiation field exactly, how- 
ever we assume that the source function S v (r) is a known 



2.2 Ray Equations 

In spherical symmetry the radiation field may be described 
as a function of radius r and the angle 6 to the local radius 
vector. A pair of variables more appropriate for describing 
the properties of the radiation field along a particular ray is 
given by the impact parameter, p = rsinf? = r(l — p 2 ) 1 / 2 , 
and the distance along the ray from the point of closest 
approach to the origin, z = rcos8 = rfi where fi = cos#. 
Radiative transfer along a ray of a given impact parameter 
p is described by the equations 

, 21 dl f{ P ,z) 

dv (4) 



dz 



■ flvir) - Xv{r)I v (p,z) 



where I^(p, z) is the specific intensity of radiation a distance 
z along the ray, in the direction of increasing z for ij" or 
decreasing z for I~ . In spherical symmetry the intensities are 
equivalently written as I„(p,z) = I v {r,fi) and I^(p,z) — 
I v {r,—fi) where fj, > 0. It is convenient to work with the 
following linear combinations of the specific intensities: 



U V {P;Z) = ~[IZ(P,Z)+I V (P,Z)], 



v v (jp,z) 



5 [tf(p>*) 



h (P,z)]. 



(5) 
(6) 



We may obtain equations for the development of u v and v v 
along a ray by summing and differencing, respectively, the 
transfer equations for ij" and 1 Z ■ If we then change variables 
to the optical depth along the ray, defined by dr„ = — Xu dz, 
we arrive at 



ot„ av 
where we have defined coefficients 
a(r) 



■y„{z) 



[1 - /i 2 + /3(rV] 



(7) 
(8) 

(9) 



and S v (r) = T)u(,1t)/xv(t) is the source function of the scat- 
tering medium, which is here assumed to be a known func- 
tion. 
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We specify a maximum and minimum radius of the 
spherical system under consideration, given respectively by 
R and Rc- (The subscript 'C here denotes the 'core', the 
central object/region of the spherical system.) The bound- 
ary conditions are: 

(a) r = R, z = 2 m a X : We assume that R has been chosen 
to be sufficiently large that there is no radiation scattered 
back into the system from this radius, i.e. I v {R,[i) = for 
ft < 0, and thus along the ray I^(p, « m ax) = where z max = 
(R 2 — p 2 ) 1 / 2 . From the definitions given by equations ([5j and 
((SJl we see that v l ,(z max ) — u„(,z ma x). We arrive at a suitable 
boundary condition by substituting this into equation (0: 



~r (2max ) 



iii/( z m»)- (10) 



(b) r — Rc, z — z m i n : For the inner boundary condition we 
need to consider two distinct cases applicable depending on 
the impact parameter of the ray chosen: 
(i) p > Rc '■ We compute the radiative transfer along the ray 
which is considered to originate from z m i n = 0. Due to the 
symmetry of the problem we may impose = I~ at z = 
and therefore v v (p, jz m in) = 0. Thus equation (O becomes 



du u (z = 0) 
dr v 



= 0. 



(11) 



(ii) p < Rc'- The ray intersects the core, and the minimum 
value of z denotes the distance along the ray at which the 
intersection occurs: z min = (R c — p 2 ) 1 ^ 2 . The boundary con- 
dition is written as: 

= Vv{Zmin) -7i>(Zmin) 5 (12) 



<9t„ 



'dv 



where we assume the function v v (Rc,fi) is known. 

If the optical depth at the core radius is sufficiently high 
then the intensity will display the angular behaviour found 
in the diffusion limit, I v {Rc,y) = Jv,c + ^H u ,cH, where J„,c 
is the core mean intensity and H Vi c is the core flux. From 
the definition of v v and the equation for /„(7?c,a0 we nn< l 



v„(Rc,n) = 



(13) 



In this limit we should take H Vi c to be a solution valid in 
the diffusion limit, if such a solution is known. If the optical 
depth at the core radius is significantly less than unity then 
we may assume a free-streaming boundary condition of the 
form 



v„(Rc,n) = H^cSnitx — 1) 



(14) 



where So is the Dirac delta function and thus the radiation 
field is directed radially outward from the source (/i = 1). 
The flux and the mean intensity are equal in this limit and 
trivially related to the source spectrum L„ by H„ t c ~ Jv,c = 
L„/(47ri?c) 2 . 

The medium is subject to a radial velocity V(r) that 
satisfies V[r) > 0, dV/dr > for all r; this ensures that 
7„ > for all (p,z). From subsequent analysis of the char- 
acteristics of the system given by equations ([7]) and we 
need to enforce an initial condition at high frequency; we 
denote this frequency by ^ max . The properties of V(r) tell 
us that every point moves away from every other point and 
all radiation intercepted at a specific point from elsewhere 
in the system is redshifted. From this argument we see that 
if ^max lies sufficiently blueward of the local line profile, such 



that no line photons ever reach ^ ma x, then 
Mr,M)L =0, 



dlu(r, fi) 



= 0. 



(15) 



(16) 



Note that equation (|15[) holds for a system interacting with 
resonance-line photons only, while the condition of zero fre- 
quency gradient given by equation (|16|) also applies in the 
more general case where we allow for continuum radiation 
as the continuum varies slowly with frequency. The equa- 
tions are solved ray- by-ray , frequency- by- frequ ency using a 
matrix equation approach l|Mihalas et alj[l97a ). 



2.3 Moment Equations 

From equation {TJ we derive angular moment equations with 
respect to fi. We obtain the following for the zeroth- and 
first-order moment equations: 



1 d(r 2 H„) 
r 2 dr 



dK v 2,K V -J V 
dr r 



d{J„~K v ) + „dK v 



dv dv 



d{H„-N„) | pdN„ 
dv dv 



(17) 



(18) 



= —XvHv 



where we have made use of the definitions of the first four 
moments of the specific intensity: 

[J V ,H V ,K„,N„] = \ f 1 d M /„(M)[l,|U,//,//]. (19) 

It is necessary to utilise additional constraints relating the 
four angular moments of the radiation intensity to obtain a 
closed system of moment equations. We assume a relation- 
ship between the moments dictated by the 'variable Edding- 
ton factor' f v {r), defined as: 



Mr) 



Mr) 



(20) 



We define a further variable Eddington factor g v {r) linking 
the first and third order moments H v (r) and N v {r): 



9"{r) 



N v (r) 
H v {rY 



(21) 



The values of fv(r) and g v (r) are dependent on the angular 
character of the radiation field at (r, v), and are determined 
by the solution of the ray equations. 

We estimate the nth-order angular moments of the in- 
tensity I v (r, /x) for n = 0, 1, 2, 3 which may be rewritten as 



Mr) 
HAr) 
K v {r) 
N v {r) 



i J I v {r,n)dn = J u u (r, (m) dp, (22) 



fj,I v (r,n)dfi 



fj,v v (r, n) dfi, 



(23) 



2~ / t JjZ Iv{r,v)^ = j H 2 U v (r,fJ>)dfJ,, (24) 
5/ / " 3/ "( r, ^ d ' U = / ^Mr,^)^ (25) 
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where the second equality in each case is reached by util- 
ising the definitions of u v (r, /i) and v v (r, fi) in terms of 
I v (r,fi) and I v (r,—fj,). We evaluate the integrals by a gaus- 
sian quadrature formula where the integral over fi is re- 
placed by a sum over different impact parameter values. For 
numerical accuracy, we split the integral into the domains 
< y, < nc and (ic < H < 1, where \i c = (r 2 + i?c) 1/2 /r 
corresponds to the core radius opening angle subtended at 
r. We then obtain f v (r) and g v (r) by evaluating the relevant 
ratios as described in equations (|20[) and (|21[) . 

We def ine the 'sphericality factor' q v (r) by the relation 
l|Auerlll97ll ) 



dln(rV) 3/„ 



dr 



fur 



(26) 



We calculate q v (r) given the variable Eddington factor f v (r) 
and the arbitrary normalisation q v (Rc) = 1 from 



q»{r) 



Rc 



exp 



r 3/„(r') 



r'U{r>) 



■ dr' 



(27) 



This quantity allows us to make a change of variable in the 
radial derivative to the dimensionless quantity X v (r) defined 
through 



dX„(r) = -Xv{r)q u (r)&r. 



(28) 



The terms on the LHS of equations (|17|l and (|18|l may then 
be expressed as 



d(r 2 H„) 
dr 



d {r 2 H v ) 
dX v 



(29) 



djfvJv) 3fu - 1 T _ -2 d{r 2 qvfvJy) 

T Jv — I Xm av ■ \ OKJ > 



Or 



dX v 



These relations and the definitions of /„ and g v may be used 
to recast the zeroth and first order moment equations as 



d(r 2 H v ) 

OXy 



8{l-f v )r 2 J v d{f v r 2 J v 



dv 

r 2 {J„-S v ), 



d(f„q v r 2 J v ) 
dX v 



■r„ 



0(1 - g v )r 2 R v _ + d(g v r 2 H„) 



dv 



dv 



(31) 



(32) 



where F v (r) = a{r)/xv[r) and once again S v (r) — 
T] v (r) /Xv(r) is the source function, although in solving the 
moment equations we account for the dependence of the 
source function on the radiation field and will not assume 
a given estimate as for the ray equations. Assuming the 
isotropic forms of the opacity and emissivity for resonance- 
line scattering, a good approximation in the comoving frame, 
the source function is given by 



S v (r) 



R(v',v,r)J v ,(r)dv' (33) 



Xu(r) <p(v,r) 

where we have allowed for the redistribution function and 
absorption profile to vary with radius due to any variation 
in the temperature profile of the medium, T(r). The integral 
is expanded as a numerical quadrature sum as described in 
Appendix lAl 

We note that boundary conditions appropriate for the 



solution of the system of moment equations are given by tak- 
ing equation (I32|l and (i) directly specifying H v (Rc) = H v ,c 
for the inner boundary, and (ii) at the outer boundary, 
specifying the ratio h v = H V (R) / J V (R) as calculated from 
the formal ray solution similarly to the variable Edding- 
ton factor computation. The equatio ns are solved using t he 
Feautrier matrix equation approach l|Mihalas et al.|[l976f ). 



3 NUMERICAL RADIATIVE TRANSFER IN 
SPHERICAL SYMMETRY II: MONTE 
CARLO 

We implement a Monte Carlo method to obtain the mean 
intensity of Lya photons and the corresponding scattering 
rate in a spherically symmetric medium with a general ra- 
dial velocity profile. It differs from existing schemes in a few 
aspects, leading to an order of magnitude increase in speed 
for a given accuracy, so we describe it in some detail. In ad- 
dition to the assumption of spherical symmetry, we sample 
the redistribution function directly and compute the energy 
density based on the path lengths traversed by the photons. 

Here we summarise the method developed and refer the 
reader to Appendix [B] for details. We impose a grid of ND 
radius values {run, Tkd-i, ri} where tnd is the inner- 
most radius at which we compute the intensity and n is the 
radius of the outer boundary. We assume the neutral hydro- 
gen is confined to r < n and thus there is no scattering for 
r > n. The algorithm, in outline, is as follows: 

(i) A photon packet is emitted at some radius r em and co- 
moving frequency i a m determined from the properties of the 
source. The packet is assigned a direction typically taken 
to be isotropic, i.e. fj, em — 2R — 1, and an optical depth 
r = — In R, where R is a random deviate uniformly dis- 
tributed ov6r ^ R 1 (drawn separately for /i e m and r), 
after which it will be scattered by the hydrogen. 

(ii) We follow the potential path of the photon packet de- 
fined by r em and /i e m, originating at A = 0, where A describes 
the path length, to the first intersection with a shell bound- 
ary and compute the distance A s and optical depth r s from 
equations (|F3 If) and (|B5|1 . respectively. 

(iii) If r > r s then the photon packet will cross the bound- 
ary, if not we skip to (iv). We update the position by set- 
ting (A + A s ) — > A and compute the comoving frequency 
between the boundaries using equations (|F32|) and (|B3[) . We 
bin the frequency and add the appropriate <5A(s) to the sum 
in equation (|B7[) for the corresponding frequency bin(s) and 
volume cell. We take (r — r s ) — > r and determine values of A s 
and t s from the current position to the next shell boundary 
along the path; this step is repeated until r s > r indicat- 
ing the packet scatters before reaching the next boundary. 
If the packet crosses the boundary at r\ then it escapes the 
H i medium and cannot be scattered back to r < n . In this 
case we return to step (i) for a new photon packet if re- 
quired. 

(iv) We assume the packet travels a further distance AA = 
(r/r s )A a along the path and add <5A = AA to the sum to 
compute J v (r) for the appropriate volume cell, where again 
the comoving frequency is binned along this section of path 
and AA must be broken into multiple values of SX if the 
range in comoving frequency is not contained in a single fre- 
quency bin. We obtain the total distance travelled by the 
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photon packet from r em before scattering from A + AA — > A. 
The radius r is updated according to equation (|B6[) with A 
in place of A max . The comoving frequency from which the 
packet is scattered, x', is determined from equation (|B3[) 
with n(r) = (fi em r em + A)/r. 

(v) The photon packet is re-emitted at r em = r with a co- 
moving frequency x em = x em (x',T) determined from the 
appropriate redistribution function, e.g. type II redistribu- 
tion (RII) with or without recoil, or even coherent scattering 
(i cm = x'). We assign an isotropically distributed direction 
ftem = 2_R — 1 and optical depth to next scatter r = — In R. 
We return to step (ii). 

This sequence of steps is repeated for as many pho- 
ton packets as are necessary to reduce the statistical noise 
in J v (r) to acceptable levels. The normalisation in equa- 
tion (|B7|) is enforced by relating the At for the Monte Carlo 
simulation to the number of photons followed, using the as- 
sumed physical properties of the source as described later 
for the particular examples we study. 



the radial derivative, is 



4 TEST PROBLEMS FOR A HOMOGENEOUS 
MEDIUM 

4.1 Static Sphere 

Diiks tra et al.l (|2006l ) obtained an analytic solution for the 
mean intensity of Lya radiation in a static uniform sphere 
allowing for RII frequency redistribut i on an alogous to that 
for a plane-parallel slab (iHarringtonl fl973) . They consid- 
ered a system of radius ii and line-centre opacity Ko, with 
a source term of unit strength, frequency distribution <j>(x) 
and arbitrary radial distribution j(r), and assumed an outer 
boundary condition applicable in the Eddington limit: 



Hu{r)\ R = -Ur)\ R . 



(34) 



They found the solution 



J(x,r) = 



6-7TKO 

16tv 2 R 



sin (A„r) 
A„r 



■ exp 




where the factors {Q n } are related to the radial distribution 
function j(r) by 



' , 2 sin(A„r) 
4Tvr -j(r)dr, 



(36) 



and the radial 'eigenvalues' {A n } are constrained by the 
outer boundary condition to the frequency-dependent val- 
ues: 



A» (a:) w — 



1 - 



l + 1.5n 1 / 2 n R(l}(x) 



(37) 



In the diffusion approximation assumed by iDiikstra et al.l 
(2006), the radiation flux is related to the average inten- 
sity by the first order angular moment equation, H(r,x) = 
— l/(3xu)dJ(r, x)/dr where the opacity is \v — \/ttko4>(x). 
The corresponding solution for the flux, derived by taking 



H(x,r) = 



A8-K 2 Rd>(x) ^ A„ 

n— 1 



Qn 



i(A n r) A n cos(A n r) 



x exp 



/ 27T A n . 3 , 

27 aK ' ' 



(38) 



IDiikstra et al.l ()2006l ) suggest a source term describing 
a thin shell at some source radius rs with the radial distri- 
bution and corresponding Q n values 



fe( 



d(t — rs 



1 sin(A n r s ) 



/7TK0 



47J-J.2 



/7TK0 



rs 



(39) 

which corresponds to a point source in the limit rs — > or 
r ^> rs. The source term is S(x,r) = Xvj{ r ) I '(4ir). The nor- 
malisation ensures the condition J dx J dV § dflS(x,r) 
I is met. 

We now consider the application of the moment and ray 
equation solution methods to this particular problem. The 
inner boundary condition is applied at a 'core radius' Rc- 
As the moment equations that we solve do not account for 
a source term, it is important that the source is contained 
within the core, Tic > rs, such that the properties of the 
source are described by the core boundary condition. This 
boundary condition is given by assuming the diffusion limit 
v v (Rc,n) = ?>\{i\H v (Rc) in solving the ray equations. In 
both cases we take the analytic solution equation (|38p for 
H u (Rc) with the Q n coefficients given by equation (|39[) . In 
estimating the source function assumed when solving the 
ray equations, it is sufficient to evaluate the integral equa- 
tion (|33|l using the analytic solution equation (|35[l as an 
estimate of the mean intensity, while the source function in 
the moment equations follows from equations ()A4|) and ()A6|I 
with e = for the coefficients. We take a linear grid in ra- 
dius spanning [Rc , R] and a linear grid of frequency values 
{xk} spanning the range over which the analytic solution is 
non-negligible. As a first step we attempt to solve the prob- 
lem in the Eddington approximation by solving the moment 
equations with f v {r) = 1/3. 

Examples of the resulting moment equation solution are 
shown in Fig. [T] We find our results agree with the analytic 
solution of Dijkstra et al. apart from differences across the 
line centre in the lower optical depth case; the moment equa- 
tion solution is the more accurate of the two as the analytic 
solution is an exact solution in the Eddington approximation 
only in the limit {ar ) l/s > 1. 

A numerical solution without assuming the Eddington 
approximation may be constructed by solving the moment 
equations with Eddington factors determined from solving 
the ray equations as described in Section 12.11 In Fig. [2] we 
compare the ray /moment equations solutions with the corre- 
sponding Monte Carlo solutions. The Monte Carlo solution 
in the upper panel of Fig. [2] was obtained according to the 
summed path-length method of Section [3] while in the lower 
panel we obtained the mean intensity at the outer boundary 
surface using the definition of specific intensity. The Monte 
Carlo solutions in each case were normalised by computing 
the simulation timescale At as required by equation (|B7|) or 
equation (|B8|I as At = N/(l photon s _1 ), which follows from 
the unit source strength imposition. The solutions were ob- 



© 0000 RAS, MNRAS 000, 000-000 



The scattering of Lya radiation in the intergalactic medium: numerical methods and solutions 7 





-20 -15 -10 



10 15 20 




-200 -150 -100 -50 



100 150 200 



-200 -150 -100 -50 50 100 150 200 



Figure 1. The mean intensity J(r,x) for Lya RII scattering in 
a static, uniform sphere in the Eddington approximation. Solu- 
tions assume a source radial dependence of the form equation d 39 D 
for rg = 5, a temperature T = 10 K, and a line-centre opacity 
kq = 100 (upper panel) and kq = 1.2 X 10 5 (lower panel). Fre- 
quency profiles are given from r = 500 (highest values of J) to 
T = 1000 (lowest values of J) in steps of 50. In each panel we 
display the solutions of the moment equations in the Eddington 
approximation (solid line) together with the cor responding an- 
alytic solution in the diffusion approximation of iDiikstra et aU 
(2006) (dashed line). 



tained by following TV = 2 x 10 5 photon packets for the case 
shown in the upper panel, and N — 2 x 10 4 photon packets 
for the lower. The frequency redistribution was treated using 
the lookup table method in the low optical depth case, while 
in the optically thick case we used the scattering atom ve- 
locity method with a bias to skip over core scatterings. From 
Fig.[2l we find that our ray/moment equation solutions agree 
closely with the Monte Carlo solutions for both low and high 
optical depths, with only small deviations from the analytic 
solution in both cases. We found no significant alterations to 
the ray/moment solutions upon re-solving the ray and mo- 
ment equations with an estimated source function obtained 
from the moment solution, suggesting that for this problem 
the combined method converges in a single iteration. 

4.2 Homogeneous Expanding Medium: 
Lya Source 

As a simple test of the ray /moment equation and Monte 
Carlo methods for a non-static, spherically symmetric 
medium, we use the problem of a L ya point source in a ho - 
mogeneous expanding H i medium (|Loeb fc RvbickiH l999L 
Loeb & Rybicki added an extra source term S(v, r) to the 
RHS of equation Q to model a point source of Lya photons 



Figure 2. The mean intensity J(r,x) for Lya RII scattering in 
a static, uniform sphere. Solutions assume a source radial de- 
pendence of the form equation d 39 1 ) for rg = 5, a temperature 
T = 10 K, and a line-centre opacity kq = 100 (upper panel) and 
hq = 1.2 X 10 5 (lower panel). In the upper panel, frequency pro- 
files are given from r = 600 (highest values of J) to r = 1000 
(lowest values of J) in steps of 50 while we display only J(R, x) 
in the lower panel. We compare our ray/moment equation so- 
lutions (solid line) with the corresponding Monte Carlo solution 
(dashed line) and the analytic solution in the diffusion approxima- 
tion (dotted line). While noisy near line centre, the Monte Carlo 
solutions agree well with the ray/moment equations solutions in 
the wings. 



at the origin: 



S(v,r) = N a 8{v-v a ) 



(47T9-)' 



(40) 



where N a is the rate of emission of Lya photons with fre- 
quency u a . They introduced a characteristic frequency v, , 
defined such that photons in the red Lorentz wing with 
Xv = tj(u — u a )~ 2 , where r\ = nHoT Q /(47r 2 ), accumulate 
an optical depth of unity in redshifting from the source with 
frequency v = v* to an observer at infinity: 



1 

a 



aY a \ a nn{z) oT Q A Q nH(0) 



5.52 x 10 



4n 2 H(z) 



(l + z) 3/2 Hz 



1/2 



(1 + Z 



,3/2 



(41) 



where we have assumed a neutral medium with tih(z) = 
7i H (0)(l + zf, H(z) ~ H nU 2 (l + z) 3/2 , as the universe 
is matter-dominated at the redshifts of interest, and a = 
H(z)v a /c. Loeb & Rybicki defined the characteristic radius 
r* as the radius at which a photon free streaming from the 
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source redshifts from v a to (u a — u t ), i.e. ar, = v t : 

v*_ _ aT a A 2 n H (.z) ^ aT a A 2 7i H (0) 
r * a ~ 4ir 2 H 2 (z) ~ 47r 2 n mJ ff2 

* 6.72(^)m P c. (42) 

Dimensionless frequency and radius variables are then given 
by v — (y a — u)/u„ and r = r/r*. In these units, a photon 
emitted at frequency £ cm will redshift over a distance r to 
u(f) — u cm + f. It is also useful to define a dimensionless 
radiation intensity / = I v /I* or J = J v /l\, where l\ = 
N a /(r*i/ t )- Loeb & Rybicki derived an analytic solution for 
the mean intensity arising from a Lya point source in the 
homogeneous expanding medium, applicable in the diffusion 
limit: 

^^iH^r^-s)- (43) 

The corresponding scattering rate per atom is 

P a = 4tr7 J J„{r)<p(v) Av = ol\P a , (44) 
where the dimensionless scattering rate P a is given by 

Here 7 = H(z)/[aXc,nn(z)] is the Sobolev parameter for a 
uniform velocity gradient H(z). The flux in the diffusion 
limit is straightforwardly derived from J: 

H(?,*) = ~¥ = — fAfW-^l (46) 

where x = r *Xf — £ -2 . Details of the derivation are pro- 
vided in Appendix [C] 

In this section, we attempt to reproduce the analytic 
solution using the moment equation solution method, be- 
fore solving the problem more generally (outwith the dif- 
fusion limit) using the ray and moment equation methods 
and comparing with the Monte Carlo solution. We continue 
to assume coherent scattering. We take our grid in radius 
to be spaced logarithmically from Rc to R. Anticipating 
results comparable to those shown by Loeb & Rybicki, we 
adopt R c = 10~ 3 and R = 10 2 . The inner boundary con- 
dition for the moment equations is given by directly speci- 
fying the flux at the boundary according to equation (146[) , 
while in the ray equations we assume the same form for the 
flux with an angular dependence given by the diffusion limit 
as in equation (|13[) . The outer boundary condition results 
from assuming no photons are scattered back into the sys- 
tem from f > R. This is technically not guaranteed to apply 
in an infinite scattering medium, but we expect no signifi- 
cant effect on the solution at distances well within the outer 
surface. Our frequency grid is subject to the restriction that 
the highest frequency (and therefore smallest value of v) lies 
bluewards of any frequency at which there are photons, in 
order to ensure compatibility with the frequency initial con- 
dition equation (|15[) . This may be ensured by taking v\ to 
lie blueward of the frequency of the least-redshifted photons 
emitted by the source as arises from free streaming, which 
requires that v\ < Rc- We use a logarithmically spaced grid 
spanning the range from v\ = to z>nf = 10 . 




r 



Figure 3. The mean intensity J(f , £>) (upper panel) and flux 
H (f , v) (lower panel) for a point emission line source in a uni- 
formly expanding homogeneous medium, from the moment equa- 
tions in the diffusion approximation (solid lines) together with 
the analytic diffusion solutions (dashed lines). Coherent scatter- 
ing is assumed. Profiles are given for log 1() v = —1.5, —1.0, —0.5, 
0.0, 0.5, 1.0, 1.5 in order of decreasing values of J or H . 



A zero-temperature medium is assumed that provides 
no Doppler broadening of the Lya line, giving a Lorentz 
profile for the opacity which in the wings takes the form 
X = 1/u 2 . The frequency redistribution function results in 
coherent scattering, so that the source function is simply 
S u (r) — J v (r). From equation (|A4|I . the redistribution coef- 
ficients for coherent scattering simplify to lZk',k.d — >• &k',k- 
When estimating the source function in order to solve the 
ray equations it is sufficient to take Sk,d ~ Jk,d and assume 
the analytic solution equation (|43l) . The analytic expressions 
equations (|43[) and (|46[) apply in the diffusion approxima- 
tion where the Eddington approximation is used, and the 
frequency derivative of the flux is assumed to be negligible. 
The former point requires f v (r) = 1/3 (note that as the ve- 
locity law is linear, /3 = 1 and the equations are independent 
of g„(r)), while the latter is enforced explicitly by setting to 
zero corresponding elements in t he matri c es use d to solve 
the problem (details are given in iHiggins! (|2012M . The re- 
sulting solutions of the moment equations for J and H are 
given in Fig. [3] They agree very well with the analytic so- 
lutions, deviating only beyond the cutoff radius where the 
finiteness of the frequency grid restricts the solution from 
following the steep decline with radius. 

A more comprehensive illustration of the utility of the 
combined ray and moment equation solution method is pro- 
vided by seeking an exact solution to the homogeneous ex- 
panding medium problem and comparing our results with 
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Figure 4. The mean intensity J(f , v) (upper panel) and flux 
H(f, v) (lower panel) for a point emission line source in a uni- 
formly expanding homogeneous medium, obtained from two it- 
erations of the combined ray /moment method with a diffusion 
solution inner boundary condition (solid lines), compared with 
the respective Monte Carlo solutions (various points). Coherent 
scattering is assumed. Profiles are given for log 10 v = —1.5, —1.0, 
—0.5, 0.0, 0.5, 1.0, 1.5 in order of decreasing values of J or H . 



the Monte Carlo solutions. We take as an inner boundary 
condition the analytic solution equation (|46[) ; this is phys- 
ically reasonable as the moment equations are expected to 
converge to the diffusion limit at great optical depths. The 
corresponding Monte Carlo solution for the mean inten- 
sity J is obtained from the summed path-length method 
of Section [3] while the solution for the flux H is obtained 
from the definition of specific intensity and the appropri- 
ate angular integration. In Fig. [4] we compare the results 
of the combined ray /moment equation solution method and 
the Monte Carlo solution method. The Monte Carlo solu- 
tions are normalised by expressing the timescale At in equa- 
tion (|B7|) and equation (|B8|) as At = N/N a . The solutions 
are based on N = 10 7 photon packets. Both the mean in- 
tensity and the flux obtained from the ray /moment method 
closely match the corresponding solutions determined from 
the Monte Carlo method. The cusps at large f for the higher 
values of log 10 v agree with those found by Loeb & Ry- 
bicki, who attributed them to causal limitation at the sur- 
face f = v. The solutions were obtained from two iterations: 
the moment solution, determined using Eddington factors 
from the ray solution, determines the source function esti- 
mate for a second ray solution, which determines the Ed- 
dington factors for a second moment solution. Upon further 
iterations we found that both methods had converged upon 
the same solution, independent of the initial estimate of the 



source function used in the first ray solution, as we tested 
by using a free streaming source function estimate of the 
form J = Sb(u — f)/(47rf) 2 . We also solved the problem us- 
ing an inner free streaming boundary condition of the form 
H = 5t>{v - R c )/(4-kRc) 2 - Except on the inner boundary, 
the solutions are identical to those produced by the diffusion 
boundary condition. 



4.3 Homogeneous Expanding Medium: Analytic 
Results for a Continuum Source 

We consider the problem of Lya scattering about a contin- 
uum point source of UV photons in the neutral hydrogen 
intergalactic medium at high redshift, prior to the large- 
scale onset of reionisation. It is assumed that the scattering 
medium may be treated as spherically symmetric, at least in 
the vicinity of t he source. The problem is then identical to 
that treated bv lLoeb fc Rvbic ki (1999), with the exception 
of the frequency dependence of the source, which is here 
taken to be spectrally flat. The source term analogous to 
equation (|40|) is 



S(u,r) = N„ 



goo 

(47tr) 2 



(47) 



where the constant N v is the rate of emission of photons per 
unit frequency. The parametrization of the radiative transfer 
equation proceeds identically to the case of a Lya source, 
with the exception of the characteristic intensity which is 
defined as = N v /r1. It is useful to generalise the frequency 
dependence of the source by assuming a limiting frequency 
Vm bluewards of the resonance frequency beyond which no 
photons are emitted: 



S(v,r) 



(47T7-) 2 



(48) 



where the step function Q(y m — v) is equal to unity for 
v < v m (or equivalently v > £> m ) and is zero otherwise. 
We show in Appendix [C] that the resulting radiation field in 
the diffusion limit is described by the solutions 



J(f,i>) 



16 



1/3 



1 



(4 7r )5/2 f7/3 



du [t + 



-2/3 



-1/2 - 

u ' e 



(49) 



3 4 



1/3 



(47r) 5 / 2 f 10 / 3 



du It + 



-2/3 



1/2 - 

u e 



(50) 



where t = — 4£ 3 /9f 2 . Note that in the case of a spectrally 
flat source with no blue cutoff frequency, v m — > — oo and the 
lower limit in each integral reduces to zero. 

The solutions are shown as a series of frequency profiles 
in Fig. [S] (Some small numerical artefacts resulting from 
the numerical integration have been interpolated over for 
the purpose of presentation.) We find a peak in both the 
mean intensity and the flux that drifts redwards further from 
the source, corresponding to the contribution from photons 
emitted by the source at the line centre. This is seen from 
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Figure 5. Analytic solution for J (upper panel) and H (lower 
panel) for a spectrally fiat source in a uniformly expanding homo- 
geneous medium in the zero-temperature diffusion approximation 
(solid lines), together with the corresponding analytic solution of 
Loeb & Rybicki for a monochromatic Lya source (dashed lines). 
Values of r are given by 10" 4 ' 5 , 10" 4 ' 2 , 10" 3 ' 9 , 10" 3 ' 6 , 10" 3 ' 3 , 
10 -30 , 10 -2,7 in order of decreasing J or H. 



a comparison of the peak frequency with that found from 
the analytic solutions for J and H for a monochromatic 
Lya source, given by equations ()43|) and (|46|) . By solving 
dJ jdv — for the analytic solution of Loeb & Rybicki, we 
find the frequency of the peak satisfies v = (3/2) 1,/3 f 2/ ' 3 . 

The most physically relevant aspect of the Lya solution 
to the 21cm signature is the Lya scattering rate. As shown in 
the upper panel of Fig. O J(x) is approximately flat across 
the line centre at radii exceeding some minimum distance 
from the source. The scattering rate then varies as the av- 
erage intensity at line centre, 



4rrcr / J u (r)y{v)dv ~ 4naI*J(r,0) (51) 



where 
J(f,0) = 



2 7/3 



3V3 (4^)5/2^:7/3 

2 7/3 
31/3( 47r )5/2jx7/3 



9f2 



t 1/6 e~ 



dt 



I I I, 



9? 2 



6' 4|£> n 



(52) 



and F(q,y) = P 



' du is the upper incomplete 



Gamma function. The solution saturates with F(7/6, y) — > 
r(7/6) ~ 0.93 for values of the argument y < 10 -2 or, for 
our problem, r < 0.1 X (2/3)j£ ; m j 3/2 . In this limit, the di- 
mensionless scattering rate P a — P a /(al^) becomes 

gdifl _ 2 7 / 3 r(7/6) 




Figure 6. Scaled scattering rate for diffusion solution with flat 
spectrum extending infinitely bluewards, P^ lfI (solid line) and 
expected result in the corresponding free-streaming limit, P^ rcc 
(dashed line). 



This form for the scattering rate applies for all r if |^ m | — > 00 
and thus applies in the absence of any cutoff in the source 
spectrum. For a finite |z/ m | the effect of the cutoff is to sup- 
press J(f, 0) for sufficiently large f: physically this corre- 
sponds to the redshifting past line centre of all photons at 
these radii. 

A useful comparison may be made with the expected 
scattering rate in the free streaming limit. The solution is 
J = Q{y — v m — f)/(47rf) 2 and the corresponding scattering 
rate for sufficiently negative v m is simply 



Anr 2 



(54) 



We show the radial dependence of the Lya scattering rate 
induced around the source in Fig. [6] The diffusion limit 
scattering rate exceeds the corresponding result in the free- 
streaming limit close to the source and decreases to below 
the free-streaming rate at f > 2 7 (r(7/6)) 3 /(3(47r) 3/2 ) ~ 
0.75. 

The size of the Lya scattering region around the source 
is given by the maximum radius at which the scattering rate 
is non-zero. This radius will depend upon v m , which we take 
to be the bluemost frequency at which the source emits pho- 
tons that are capable of redshifting into the Lya line. Typi- 
cally this cutoff will be the Ly/3 frequency, Vp = (32/27)^ a , 
as photons emitted bluewards of Ly/3 will redshift into the 
Ly/3 or higher-energy resonance lines and do not produce 
Lya photons (except as products of radiative cascades which 
we do not consider here). Assuming a Lya source term of the 
form of equation (I48p . the scattering rate found in the diffu- 
sion limit given by equation (|52l) is found to decrease below 
one-tenth of its value in the absence of a cutoff when the sec- 
ond argument y of the upper incomplete Gamma function 
satisfies y > f ~ 2.5, or 



f >(2/3)/ 1/2 |^| 3/2 ~1.05|i> /3 | 3/2 - 



(55) 



31/3 (47r)3/2 fr/3- 



(53) 



Outside of the diffusion regime, an upper limit for the size 
of the Lya scattering region is given by free streaming: all 
photons will have redshifted past the line centre at radii 

f > \v (j \. 

The number of scatters a Lya photon undergoes before 
escaping the IGM may be computed as the ratio of the total 
rate of scatters through a region of radius r max and the 



© 0000 RAS, MNRAS 000, 000-000 



The scattering of Lya radiation in the intergalactic medium: numerical methods and solutions 11 



emission rate: 

iVscatt = iVT^TT 



drr n H (z)P a (r) 



12 2/3 (1_ 
(4tt)V2 ^6 

-1 'max 



— 1 — l 'max 
' hor ^"hor "I - ~ 

' hor 



''ho 



(56) 



7 



where equation (|53|) is used for < f < 1 and equation (|54[) 
for f > 1. Here rh or = (5/27)c/H(z) is the 'horizon' dis- 
tance a photon emitted just longward of the Ly/3 resonance 
frequency may travel before redshifting into the local 
Lya resonance frequency, and the total Lya photon pro- 
duction rate is taken as N a = N v (vp — v a ) = (5/27)N v v a . 

For the Wouthuysen-Field effect to compete with the 
CMB and couple the spin temperature to the gas kinetic 
temperature, a critical thermalization Lya scattering rate 
of Pth = 27^1ioTcmb/4T, is required, where Aw ~ 2.85 x 
10 _15 s _1 is the 21-cm transition rate, T» = hvio/ks — 
0.068 K and Tqmb is the temperature of the Cosmic Mi- 
crowave Background (|Madau et al.l ll997T ). For a contin- 
uum source, from equations (|5ip and (|53|l the critical ther- 
malization photon production rate between the Lya and 
Ly/3 resonance line frequencies is 



4.39 x 10° 



1 + z 
11 



7/3 - 
r Mpc S 



(57) 



where vm P c is the distance from the source in megaparsecs. 
Here we adopted properties corresponding to the mean in- 
tergalactic medium at z = 10, taking the baryon and dark 
matter density parameters flth 2 = 0.022 and Q, m — 0.27, 
respectively, and Hubble con stant Hp = 70kms~ 1 Mpc" 1 , 
consistent with the results of iKomatsu et aD (|201lh . These 
values give u* = 1.25 x 10 13 [(1 + z)/ll] 3/2 Hz and r* = 
1.12 Mpc (independent of redshift). The Sobolev parameter 
is 7 = v a H/{acnw) ~ 1.27 x 10" 6 [(1 + z)/ll]~ 3/2 . Thus 
a source with a continuum luminosity of ~ 2 x 10 11 Lq 
would be able to couple the spin temperature to the gas 
kinetic temperature out to a distance of 1 Mpc. This is far 
less demanding than for an emission-line source. From equa- 
tions 14411 and (I45[) . the required Lya photon production rate 



N" 



1.38 x 10 b 



1 + z 
11 



11/3 -1 
r Mpc S 



(58) 



An emission-line source with a Lya luminosity as great as 
10 12 Lq would be able to couple the spin temperature to the 
gas temperature only out to a distance of ~ 90 kpc. Indeed, 
at this distance, if the source were a radio-loud AGN, its 
radio emission wo uld likely dominat e the coupling of the 
spin temperature l|Madau et al.ll 19971 ') . 



4.4 Homogeneous Expanding Medium: Numerical 
Results for a Continuum Source 

Outside of the zero-temperature regime treated by Loeb & 
Rybicki for a monochromatic source, it is necessary to spec- 
ify the physical state of the scattering medium. We adopt 
properties corresponding to the mean intergalactic medium 
at z = 10 as above, with an assumed gas temperature 
T = 10K. 

In solving the moment equations, the initial condition in 
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Figure 7. Diffusion approximation moment soiution for J in 
a uniformiy expanding homogeneous medium assuming coher- 
ent scattering and a Voigt line profiie (solid line), together with 
the corresponding analytic solution in the wing approximation 
(dashed line). The frequency dependence is shown as a function 
of x assuming a temperature T = 10 K for values of f given by 

io- 4 5 , id- 4 - 2 io- 3 - 9 , io- 3 - 6 , icr 3 - 3 10- 30 . io- 2 - 7 , io- 2 - 4 . 



frequency is enforced by taking our frequency grid to start 
from v\ < £ m , i.e. the frequency grid must extend blue- 
wards of the source emission cutoff frequency. Unless other- 
wise stated, we have taken s m = — (i/»/Ai^>)i' m = 1000 for 
all solutions. We use two frequency grids, one with a small 
frequency spacing across the line centre for x lying in the 
range [— x max , a: max ] where a; ma x = 60 (unless stated oth- 
erwise), and a second grid with a larger frequency spacing 
across [x max ,xi]. We adopt the analytic solutions to satisfy 
the inner and outer boundary conditions for the moment 
equations in the diffusion approximation. We solve the mo- 
ment equations with the addition of the full Voigt profile 
form for \.v with an opacity given by x = V^r r Ko<f>v(x). 
The diffusion limit is enforced by setting f v (r) = 1/3, and 
the appropriate matrix elements representing the frequency 
derivative of the flux in the solution method are set to zero 
to match the diffusion equation. Our moment solution for J 
is shown in Fig. [7] It nearly coincides with the analytic so- 
lution, demonstrating that the assumption of the wing form 
of the opacity does not affect the solution in this approxi- 
mation. A similarly identical solution was found for the flux. 
As a consequence, the solutions as a function of frequency 
are very insensitive to the temperature of the medium in the 
coherent scattering approximation. 

We solved the problem using the ray/moment method 
applied to the exact forms of the ray and moment equations 
for the uniformly expanding homogeneous medium problem. 
Again, we used the Voigt profile form for x an d assumed co- 
herent scattering. We assumed the diffusion limit solutions 
given by equations (|49[) and (|50[) in specifying the inner 
boundary condition and the initial estimate of the source 
function used in solving the ray equations. Note that for 
this and subsequent problems where the conditions do not 
match those that determine the solution used for the bound- 
ary condition, the ray /moment solution at the boundary is 
subject to an unphysical constraint and we should not as- 
sume it is accurate on the boundary. 

Our ray/moment solution is displayed in Fig. [S] For 
distances sufficiently far from the source we find the 
ray/moment solution is a close match to the diffusion ap- 
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Figure 8. Coherent scattering solution for J obtained from 
ray and moment equations (solid lines), compared with the an- 
alytic solutions in the zero-temperature diffusion approxima- 
tion (dashed lines). Frequency profiles correspond to radii f = 

icr 4 - 2 , ur 3 - 9 , icr 3 - 6 , icr 3 - 3 , 10- 30 , icr 2 - 7 , nr 2 - 4 . 
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Figure 9. RII redistribution solution for J obtained 
from ray and moment equations (solid lines), compared 
with the corresponding solution for coherent scatter- 
ing (dashed lines). Frequency profiles correspond to radii 
f = IO" 42 , IO" 39 , IO" 36 , IO" 33 , IO" 30 , io- 2 - 7 , io- 2 - 4 . 



proximation solution. Closer to the source we find both 
solutions match reasonably well across the line centre, 
while away from line centre there is some deviation. The 
ray/moment solution approaches a flat profile closer to the 
free-streaming solution. This is an expected consequence of 
the negligible optical depth out to these radii away from line 
centre. We have checked the Eddington factor and found 
f v — 1/3 except in the wings, with the deviation from 1/3 
matching the deviation in the wings between the two solu- 
tions shown in Fig. [8] We found that the flux is altered from 
the analytic solution to a similarly limited extent. 

The accuracy of the solutions may be checked by use 
of an integral constraint on the number of photons ob- 
tained from the radiative transfer equation. It may be de- 
rived from the zeroth-order moment equation obtained from 
equation (|C1[) with V = Hr and 8(9) — > 1 for a continuum 
source in a homogeneous expanding medium by multiplying 
it by f 2 , integrating in radius from f = to R and in fre- 
quency across the line centre from — £ max to -\-P maK , causing 
the terms corresponding to the opacity and emissivity to 
cancel out due to radiative equilibrium. The result is 



[j(f, 5max) - J(f, -5max)] df 



+R 2 I H(R,i>)Ai> = 2 " 



(59) 



(4tt) 2 



We have verified that all the solutions examined in this sec- 
tion satisfy this constraint. 

The assumption of coherent scattering is valid for pho- 
tons scattering at frequencies far from line centre, but as 
we examine photons redshifting across the line centre such 
an assumption is not physically motivated. We now treat 
the added effect of Doppler shifts due to the thermal- and 
recoil-induced atomic velocities. The thermal velocities of 
the atoms are responsible for the RII frequency redistribu- 
tion function for Lya scattering. For convergence with re- 
coils, we found it necessary to increase x max to 200. We first 
consider the addition of RII redistribution in isolation by ig- 
noring the effect of recoils. Ray/moment equation solutions 
are obtained similarly to the coherent solution with the RII 
form for the coefficients TZk',k,d with e = as given in equa- 
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Figure 10. RII redistribution solution with recoil obtained 
from ray and moment equations (solid lines), compared with 
the corresponding solution for RII redistribution without re- 
coil (dashed lines). Frequency profiles correspond to radii f = 
IO" 4 2 , IO" 3 - 9 , IO -3,6 , 10 -3 - 3 , IO -3 , 10 -2 - 7 , IO" 2 - 4 . 



tion (|A6I) in Appendix [X] We display our solution for J in 
Fig. [5] There is some redistribution of photons across the 
line centre relative to the red peak found in the coherent 
scattering case, resulting in a small boost in the mean in- 
tensity across the line centre over a limited range in radii. 
The difference in the spectrum from the coherent scattering 
case is not particularly significant at large radii. 

We add the effect of atomic recoil by seeking 
ray/moment solutions without setting the recoil parameter e 
to zero. The resulting solution is displayed in Fig. llQI The ex- 
pected Boltzmann distribution gradient is recovered across 
the line centre. 

We also apply the Monte Carlo method described in 
Section [3] and Appendix [B] to the problem. We assume a 
continuum source that emits N u photons per second per Hz 
across the frequency range [aimin, imai] where N v is indepen- 
dent of v, 



N v Al> = -r- X 

At 



Ax 



(60) 



where Av is the frequency bin width and the factor in brack- 
ets is the fraction of the N photons emitted within a single 
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Figure 11. Monte Carlo solution (dashed lines) for coherent 
scattering, together with the corresponding ray/moment equa- 
tion solution (solid lines). Frequency profiles correspond to radii 
f = IO" 42 , IO" 39 , IO" 36 , IO" 33 , IO" 30 , io~ 2 - 7 , io- 2 - 4 . 
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Figure 12. Scattering rates derived from numerical integration 
of ray/moment equation method (solid line) and Monte Carlo so- 
lution (long-dashed line), both for coherent scattering; normalised 
by the analytic result in the diffusion limit. 
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Figure 13. Scaled scattering rate for the IGM at z = 10 with 
an assumed temperature T = 10 K, derived from numerical inte- 
gration of the scattering solution obtained from the ray/moment 
equations (solid line) and the analytic diffusion solution (dashed 
line, nearly coincident with the solid line in the upper panel), for 
source cutoff frequency z m = 100 and coherent scattering (upper 
panel) and x m = 1.4 X 10 5 including RII redistribution and recoil 
(lower panel). For comparison we also show p^ ee (dotted line). 



frequency bin. In the case of a source cutoff at x m , we take 
i mal = x m . The solution is normalized by using At from 
equation (|60l) in equation (|B7[) . The resulting dependence 
on the source strength N v is scaled out when expressing the 
solution in terms of J. We obtained a solution for coherent 
scattering using N — 2 x 10 5 photon packets. The solution is 
displayed in Fig. 1111 It closely matches the ray/moment so- 
lution. The corresponding scattering rate is shown in Fig. 1 121 
It closely resembles the ray/moment solution scattering rate, 
which itself is nearly coincident over this range in radius with 
the diffusion approximation power law of equation (|53[) . The 
noise level of the Monte Carlo solution is a few to several 
percent per radial bin, decreasing with increasing radius. 
Two hundred logarithmically spaced radial bins were used 
to cover the range —4.5 < log 10 f < —1.5. 

We have determined the scattering rate profile from 
the ray/moment solution. Fig. 1131 shows two examples with 
different values of the source upper frequency cutoff !/ m . 
The upper panel has x m = 100 and thus |£> m | ~ 10 -1 " 6 . 
The 'diffusion-limited' effective radius of the Lya scattering 
region, as given by equation (|55p for general v m , is f = 
1.05|i> m | 3//2 — 10~ 2 ' 3 , in approximate agreement with the 
figure. 

The lower panel has x m = xp, where xp = 



5^ q /(27Ai^d) — 1.4 x 10 5 denotes the source cutoff corre- 
sponding to the Ly/3 frequency. This corresponds to —v m ~ 
10 1 ' 6 . Photons with larger values of — v would not be able 
to redshift into the Lya resonance at any greater distance, 
as they would encounter the Ly/3 or a higher order Lyman 
resonance en route. The limiting distance the Lya photons 
may travel is fhor — 38[(1 + z)/ll] -3 ^ 2 (the Lya horizon). 
For r < 1, the diffusion limit applies and the scattering rate 
takes on the approximate f~ 7 ^ 3 scaling expected. At very 
small values, f < 10 -3 , frequency redistribution substan- 
tially modifies the scattering rate, as shown in Fig. Q3] As 
discussed in Appendix[C] the diffusion approximation breaks 
down at r > 1 and the scattering rate takes on a profile closer 
to the free-streaming v alue l/47rf 2 , a s found in the Monte 
Carlo computations of ISemelin et al.l l|2007l ). The scatter- 
ing rate becomes vanishingly small at f > 29, as shown in 
Fig. [T3J This is slightly shorter than the Lya horizon dis- 
tance by the factor 0.77 and appears to be a consequence 
of the causal limitation of the radiation field (see Loeb & 
Rybicki 1999). 

The 1/r 2 behaviour of the mean intensity at large dis- 
tances has an interesting implication for the global radiation 
field produced by a uniform distribution of sources of spatial 
number density no. The total scattering rate is limited by 
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Figure 14. Upper panel: Ratio of scattering rate derived from 
numerical integration of the appropriate ray /moment solution, to 
the analytic diffusion limit result P^ lS , for RII redistribution with 
recoil (solid line), RII redistribution without recoil (dashed line) 
and coherent scattering without recoil, which coincides with the 
RII redistribution without recoil solution for f > 10 -2,5 (dotted 
line). Also shown is the free-streaming result (short-dashed line). 
All solutions assume T = 10 K and x m = 1.4 X 10 5 . Lower panel: 
Ratio of scattering rate assuming RII redistribution with recoil to 
that assuming RII redistribution alone (solid line), compared with 
the estimate S a = 0.78 (dashed line) of the suppression factor 
in a uniformly expanding homogeneous me dium for a uniform 
radiation field l IFurlanetto &: Pritchardll2006h . 

the maximum distance f max the Lya photons travel: 

P* ot = 47v{n r 3 ») / rm " dfr 2 -^— ~ (n r 3 )f max . (61) 
Jo 4ttH 

The maximum distance is given by the causal limitation 
radius f max = f causa i ~ 29 < Fhor- Given that the scatter- 
ing rate at large distances somewhat exceeds l/Aitr 2 within 
the horizon, however, taking the maximum radius to be the 
horizon radius gives a very good approximation to the total 
scattering rate, and corresponds to AT sca tt — 7~ ■ 

It is instr u ctive t o compare this result with the esti- 
mate of iFieldl (|l959al) . who expressed the scattering rate 
as the product of the production rate of Lya photons per 
neutral hydrogen atom and the number of scatters N sca , t t 
a photon undergoes before it redshifts sufficiently far from 
line centre to escape: P a — (noN a /nn)N sca ,u- He argued 
that for uniformly distributed sources in a homogeneous 
and isotropic medium, AT sca tt — 7 _1 (cf. Higgins & Meiksin 
2009). Allowing for all photons emitted between the Lya and 
Ly/3 frequency resonances, this corresponds to the dimen- 
sionless scattering rate P^ lcld = (nor 3 )fhor- For a uniform 
radiation field, the total scattering rate is thus again given 
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Figure 15. Coherent scattering solution obtained from ray and 
moment equations for the overdense shell given by equation 1021) 
(solid lines), compared with the corresponding solution for a uni- 
form density njj(^) (dashed lines). Frequency profiles correspond 
to radii f = 10" 4 2 , IO" 3 9 , 10" 3 - 6 , 10" 3 3 , 10" 3 , 10" 2 7 , 10" 2 4 . 
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Figure 16. Monte Carlo solution (dashed lines) for coherent 
scattering with the density profile specified in equation H62II , 
compared with the corresponding ray /moment equation solu- 
tion (solid lines). Frequency profiles correspond to radii r = 

io- 4 - 2 , io~ 3 - 9 , io~ 3 - 6 , io- 3 - 3 , io- 3 , io- 2 - 7 , io~ 2 - 4 . 

by equation (|61|) . with the maximum radius taking on the 
natural value f max = fhor- 

Allowing for recoils suppresses the scattering rate rel- 
ative to the RII case without recoil, as shown in the lower 
panel of Fig. 1141 The suppre ssion factor S a computed in the 
diffusion approximation by iFurlanetto fc Pritchardl (|2006l ) 
for an isotropic and homogeneous distribution of uniform 
sources in a uniformly expanding homogeneous medium is 
given by S a — 0.78. This is close to the average value of 
the radius-dependent suppression factor determined from 
our solutions outwith the diffusion approximation. 

5 APPLICATIONS TO INHOMOGENEOUS 
MEDIA 

5.1 Overdense Shell 

As an example application for an inhomogeneous density 
profile, we examine the continuum source Lya scattering 
problem for a scattering medium undergoing Hubble expan- 
sion but including an overdense shell between 10~ 3 ' 5 < r < 
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Figure 17. Scaled scattering rate derived from numerical inte- 
gration of ray/moment equation solution (solid line) and Monte 
Carlo solution (long-dashed line) for coherent scattering in a 
medium with an overdense shell given by equation Ki-'l . For com- 
parison we also show P^ l(l (short-dashed line). 
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Figure 18. Velocity profile given by the quadratic relation of 
equation (1636 with R mln = 10 -5 , -Rmax = 10 -2 and V max / m ; n = 
^^max/min (solid line), compared with the Hubble law V(r) = 
Hr (dashed line). 
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(62) 



The IGM temperature is T = 10 K. This affects the 
spatial dependence of the opacity coefficients \k,d for the 
ray/moment equations method. We show the ray/moment 
solution to the problem for coherent scattering in Fig. 1151 
while the corresponding scattering rate is denoted by the 
solid line in Fig. 1171 Compared with the spectrum for a uni- 
form density scattering medium, the frequency peaks at the 
overdense radii are displaced further redward, while the scat- 
tering rate across this range is boosted. At the outer edge of 
the shell the scattering rate dips below the uniform density 
value. 

The corresponding Monte Carlo solution to the problem 
is shown in Fig. 1161 (N — 2 x 10 5 photon packets were used.) 
Significantly more time was required to follow the photon 
packets through the overdensity compared with the homo- 
geneous density case. The corresponding scattering rate is 
displayed in Fig. [T7] 

5.2 Quadratic Velocity Profile 

In this section we consider a continuum source in a medium 
of uniform density, however we move beyond the assumption 
of a simple Hubble velocity profile V = Hr. To examine 
the specific effects of a non-linear velocity profile alone, we 
consider the case of a quadratic velocity law. We parametrize 
the radial velocity profile V(r) for i? m i n < r < i? ma x by 

V(r) = ( Knax - V- min ) ( - r - Rn ™ ) +V min . (63) 



Rn 



We show both the quadratic velocity profile given by equa- 
tion (|63|) and the corresponding Hubble profile in Fig. 1181 
where we have taken V max = V(-Rmax) = HR m!lx and 
Vmin = V(7?min) = -ff-Rmin- Using this velocity profile we ob- 
tain Lya radiative transfer solutions from the ray/moment 
method. The IGM temperature is T = 10 K. The solution as- 
suming coherent scattering is displayed in Fig. 1191 In order to 
isolate the effect of the velocity profile, the solution is com- 
pared with the equivalent ray /moment solution for Hubble 
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Figure 19. Coherent scattering solution obtained from 
ray and moment equations with the quadratic veloc- 
ity profile of equation ( 1 63 D (solid lines), compared with 
the corresponding solution for a Hubble velocity pro- 
file (dashed lines). Frequency profiles correspond to radii 
f = lO" 41 , lO" 3 - 8 , lO" 3 - 5 , lO" 3 - 2 , lO" 2 - 9 , lO" 2 - 6 , lO" 2 - 3 . 




Figure 20. Monte Carlo solution (dashed lines) for coherent 
scattering with the quadratic velocity profile of equation 1 163 II . 
compared with the corresponding ray /moment equation solu- 
tion (solid line). Frequency profiles correspond to radii f = 

io- 41 , lo- 3 - 8 , lo- 3 - 5 , lo- 3 - 2 , io- 2 - 9 , icr 2 - 6 , io- 2 - 3 . 
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Figure 21. Scaled scattering rate for a uniform medium with a 
velocity profile given by equation (163 b . derived from numerical 
integration of the ray /moment equation method (solid line) for 
coherent scattering, together with the corresponding scattering 
rate obtained from the Monte Carlo method (long-dashed line). 
For comparison we also show P^ lS (short-dashed line) given by 
equations H51H and II52H . The scattering rate decreases at large f 
near the outer boundary. 
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flow. The solution presents the same qualitative features as 
the coherent scattering solution for the Hubble velocity pro- 
file, in particular the peak in intensity that drifts redwards 
from line centre further from the source, although signif- 
icant quantitative differences are apparent. The frequency 
displacement of the peak increases more slowly with radius 
than in the Hubble velocity case, a result of a reduced ve- 
locity across the range in radius which is less efficient in 
redshifting photons. 

The computations expended by the Monte Carlo code 
are significantly reduced if we may assume the photon occu- 
pies only a single frequency bin between scattering or bound- 
ary crossing events, which we take to be that corresponding 
to the 'final' frequency prior to scattering. In a Hubble ve- 
locity field the comoving frequency changes with path length 
A as x = x arn — (H/b)\, and the approximate treatment is 
accurate in the limit that the path length between scattering 
or boundary crossing events satisfies 5X <C bAx/H, where 
Aa; is the frequency bin size, a condition which is certainly 
satisfied in our examples for distances between scattering 
events at line centre important in determining the scattering 
rate. For a general velocity profile this condition is replaced 
by 8X <C 6A:r/|dV/dr| , and so this approximate treatment 
is valid as long as V(r) does not change too rapidly. We 
obtained a Monte Carlo solution for the coherent scattering 
problem which is compared with our ray /moment solution 
in Fig. I20I The Monte Carlo scattering rate is compared with 
that determined from the equivalent ray /moment solution in 
Fig. [5T] Well interior to the outer boundary, the scattering 
rate is boosted by up to an order of magnitude compared 
with the expected rate for a Hubble velocity profile, given 
approximately by P^ lfi ' . Near the outer boundary, where the 
velocity gradient well exceeds the Hubble constant, the scat- 
tering rate dips below P^ lS . 



Figure 22. Density profile (solid line, upper panel) and velocity 
profile (solid line, lower panel) given by equations H64I I and I I65II 
respectively with amplitude Ao = 3 and compared with the cor- 
responding uniform density or Hubble velocity profile at z = 10 
(dashed lines). 



5.3 Spherical Perturbation in Density and 
Velocity 

In this section we examine the scattering rate resulting from 
a perturbation to the uniform density and Hubble velocity 
profiles of a homogeneous expanding medium. We assume a 
spherically symmetric overdensity of the form 



8(r) = Aojo(fcr) 

=> n H (r) = nn(z)[l + A j (kr)] 



(64) 



where jo{x) = (sin x) jx is the zeroth-order spherical Bessel 
function, Ao is the amplitude of the perturbation and the 
value of k is taken to be 2n/R, where R is the radius of the 
outer boundary. A restriction on the amplitude follows from 
the physical requirement that 8 > — 1, which then requires 
that Ao < 5. The corresponding self-consistent perturba- 
tion in the peculiar velocity follows from the linear regime 
equation 8 — — (1/H)X7 ■ v p which in spherical symmetry is 
solved for 8 = Aojo(fc»") and the boundary condition that 
the velocity vanishes at r = by 



k 

V(r) =Hr + v p = H{z) 



— ji(kr) 



(65) 



where ji(x) — (sin x)/x 2 — (cosx)/x is the first-order spher- 
ical Bessel function. 

We solve the radiative transfer problem for a continuum 
source subject to a scattering medium with the perturbed 
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Figure 23. Scaled scattering rate for a spherical perturbation 
with the density profile equation 1641 and velocity profile equa- 
tion 151}, with k = 2-k/R, R = 10°' 5 and x m = 6 X 10 4 , for 

Ao = 0.5 (solid line) and Ao = 2.9 (dashed line). The scattering 0.9 
rates are derived from numerical integration of the ray/moment 

equation method for coherent scattering. For comparison the rate ra g 

in the diffusion approximation for a uniformly expanding homo- 
geneous medium is also shown (short-dashed line). 
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Figure 24. Scaled scattering rate for a spherical perturbation 
with the density profile equation H64t and velocity profile equa- 
tion (65j, with k = 2n/R, R = 10 ' 5 and x m = 6 X 10 4 , as a 
function of the perturbation amplitude Ao and normalised by the 
value for Ao = 0, at different radii given by f = 10 -3 (solid line, 
P Q (0) = 7.2 X 10 5 ), 1CT 2 (long-dashed line, P„(0) = 3.5 X 10 3 ), 
lO" 1 (short-dashed line, P a (0) = 17) and 10° (dotted line, 
Pa(0) = 0.11). The scattering rates are derived from numeri- 
cal integration of the ray/moment equation method for coherent 
scattering. 

density and velocity profiles of equations (I64|l and (|65|l . re- 
spectively, using the ray/moment equation solution method 
assuming coherent scattering in a medium with temperature 
T — 10 K. We note that a further restriction on the pertur- 
bation amplitude Ao follows from requiring a monotonic ve- 
locity field, with V(r) > and V'(r) > for all r, necessary 
to apply the moment equation solution method. The result- 
ing constraint is Ao < 3, corresponding to a vanishing ve- 
locity gradient at small r. While such large values of Ao are 
no longer in the linear regime, used to derive equation (I65|) . 
the form permits an exploration of the effect of the shape of 
the velocity field on the transport of the Lya photons. 

The hydrogen density and velocity profiles correspond- 
ing to our assumed cosmological parameters at z = 10, an 
outer boundary R — 10 0,5 r* ~ 3.5 Mpc and a perturbation 



0.6 1 — ■ — — ■ — ^ — ■ — — ■ — ' — ■ 1 

1 10 100 1000 10000 

X 

Figure 25. Eddington factors /„ (upper panel) and g v (lower 
panel) for a spherical perturbation with the density profile equa- 
tion II64I I and velocity profile equation H65H for Ao = 2.9, derived 
from numerical integration of the ray/moment equations assum- 
ing coherent scattering. The factors are shown at log 10 f = 
(solid line), log 10 f = —0.5 (long-dashed line), log 10 r = —1.0 
(short-dashed line), log 10 f = —1.5 (dotted line) and log 10 f = 
—2.0 (dot-dashed line). The factors take on the values /„ ~ 1/3 
and g v > 3/5, close to the values expected for the linear approx- 
imation l v ~ J„ + 3H v /j,, across the line centre, and f v =g v = \ 
far in the wings expected for free-streaming radiation. 

with Ao = 3, are shown in Fig. [22] The scattering rate as a 
function of distance from the source is shown in Fig. 1231 for 
Ao = 0.5 and 2.9. While a small amplitude perturbation has 
a correspondingly small effect on the scattering rate, the rate 
is substantially boosted for Ao < 3. The dependence of the 
scattering rate on the perturbation amplitude for Ao < 3 is 
shown at various radii in Fig. 1241 The scattering rate scales 
linearly with the amplitude for small values, Ao < 1, but 
increases more rapidly with larger values. It grows exponen- 
tially as Ao — > 3, a result of a vanishing velocity gradient for 
r — > (cf. Figs ll9l and l21[) . Clearly a small velocity gradient 
near the source is able to substantially boost the scattering 
rate over a wide range of radii compared with the case of 
uniform Hubble flow. 

The solutions obtained from the solving the 
ray/moment equations have generally justified the lin- 
ear approximation I v ~ J„ + 3H u fi in fi, for which the first 
two Eddington factors are f v = 1/3 (the Eddington ap- 
proximation) and g v — 3/5. For a linear velocity field, only 
f v is required, as the terms involving N v in equation (JTHJ 
cancel. The perturbed flow here allows a check on g v . The 
Eddington factors are shown in Fig. 1251 They take on the 
values /„ = 1/3 and g v > 3/5 over a broad frequency 
region across the line centre, and the values f v = g v = 1 
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far in the wings, as expected for free-streaming radiation. 
On the outer boundary, h v > 0.7 over a similarly broad 
region across the line centre, allowing closure of the moment 
equations. The value corresponds to a single dominant 
escaping stream of radiation with fi ~ 2~ 1 ^ 2 . This would 
correspond to f v = g v = 1/2, close to the computed values 
on the outer boundary of /„ ~ 0.5 and <?„ > 0.6. 



6 CONCLUSIONS 

The 21 cm signatures of the first luminous objects are de- 
pendent upon the Lya scattering rate, which is necessary to 
decouple the H i spin temperature from the CMB tempera- 
ture through the Wouthuysen-Field effect in the diffuse IGM 
during the early stages of reionization and so render the hy- 
drogen detectable against the CMB via 21cm observations. 
Realistic studies of Lya scattering caused by the first sources 
during the EoR will need to make use of three-dimensional 
density and velocity fields derived from cosmological simu- 
lations, and treat the time-dependent radiative transfer of 
photons contributed from multiple sources of finite lifetime, 
the detailed properties and spatial distributions of which 
are largely uncertain and will have to be extrapolated from 
lower-redshift observations or estimated from galaxy forma- 
tion simulations. Together with the extensive effort currently 
underway in removing radio foregrounds and isolating the 
signal of the EoR, these detailed analyses will be critical 
in predicting and interpreting future 21 cm tomographic ob- 
servations made by existing and future radio interferometric 
arrays, such as LOFAR and SKA. 

We develop two methods for solving the radiative trans- 
fer equation for resonance line photons in spherically sym- 
metric systems, and apply them to idealized problems rel- 
evant to the 21cm signature of the IGM. We consider five 
classes of problems corresponding to a point source emit- 
ting either emission line or continuum radiation in an ex- 
panding medium, allowing for a uniformly expanding homo- 
geneous medium, a uniformly expanding medium with an 
overdense shell around the source, a homogeneous medium 
with a quadratic velocity profile, and a medium with a 
self-consistent density and velocity perturbation around the 
source. A static medium test problem is also treated. Since 
these problems may serve as useful tests of Monte Carlo 
schemes, we provide their solutions in machine-readable for- 
mat (Appendix [Dll. 

Following [Mihalas et~afl (|l975h and iMihalas et~ai1 
(|l976h . the first method is based on solving the ray and an- 
gular moment forms of the comoving frame radiative trans- 
fer equation in a spherically symmetric medium subject to 
a monotonic velocity profile V[r), where V > and V' > 
for all r (see Section [2. 1 [) . We described general boundary 
conditions required by these methods in order to treat the 
radiative transfer of Lya scattering in a cosmological con- 
text. We also obtained an efficient prescription for the finite 
difference representation of the source function for partial 
frequency redistribution based on the diffusion approxima- 
tion. 

We compare our results with Monte Carlo solutions to 
the radiative transfer equation using an implementation that 
determines the Lya radiation mean intensity from the ac- 
cumulated path lengths of the photons within a given vol- 



ume and frequency range. This greatly improves the accu- 
racy over a technique that uses spatial and frequency bin 
crossings to estimate the intensity. For 200 logarithmically- 
spaced grid zones covering three decades in radius, the noise 
level in the scattering rate is kept under 10 per cent per ra- 
dial bin using 2 x 10 5 photon packets. The code computes 
the Doppler redistribution of the frequencies upon scattering 
both directly and by interpolating on the RII redistribution 
function. The latter speeds the computations by a factor of 
a few, but requires a grid tailored to a specific application, 
so lacks generality. 

In Section [4] we validated our schemes using vari- 
ous test problems for resonance line photon scattering in 
spherical symmetry, includin g: (i) Lya scattering in an op- 
tically thick static sphere (|Diikstra et al.l 120061 ) and (ii) 
the Lya scattering halo of a monochromatic Lya point 
source in a uniformly expanding homogeneous medium 
|Loeb fc Rvbickil 1 19991 ). We found the ray and moment 
method results agreed with the analytic solutions, at least 
in the optically thick regime where the analytic solutions 
are valid, while in case (ii) our solution matched the Monte 
Carlo solution found by Loeb & Rybicki. A final test problem 
was described for a continuum source in a uniformly expand- 
in g homogeneous medium . We utilised the analytic solution 
of lLoeb fc Rvbickil (|l999u , valid in the zero-temperature dif- 
fusion approximation, as a Green's function to obtain the 
equivalent solution for a flat source spectrum in the diffu- 
sion approximation. The solution gives a radial profile for 
the Lya scattering rate that varies as r _7//3 , steeper than 
the free-streaming dependence r -2 , and the same depen- 
dence found by other authors in Monte Carlo studies of 
this problem. Our analytic solution also allows for a maxi- 
mum frequency cutoff from the source; both of our numerical 
methods recover the effects of the cutoff. The photon pro- 
duction rate between Lya and Ly/3 required to couple the 
spin temperature to the gas kinetic temperature is found to 
be Nl% ~ 4.39 x 10 55 [(1 + z)/ll]r 7 ^ c s"\ where r Mpc is 
the distance from the source in megaparsecs. This is less de- 
manding than for an emission-line source by several orders 
of magnitude at a distance of 1 Mpc. 

We solve the test problem of a continuum source in a 
uniformly expanding homogeneous IGM numerically, out- 
side of the diffusion limit, using our ray/moment and Monte 
Carlo methods for spherical symmetry. In suitable spher- 
ically symmetric problems such as this test problem, the 
combined ray/moment equation solution method may be 
used to quickly produce noise-free results, in contrast to 
the Monte Carlo approach for a practical number of pho- 
ton packets. For more general situations, the Monte Carlo 
method is more readily extended to treat cartesian grids 
with general configurations of sources and arbitrary density, 
temperature and peculiar velocity fields within the scatter- 
ing medium. 

For precision results, however, a grid-based scheme such 
as the method presented here may be desirable. In this case, 
the diffusion approximation may be used as an inner bound- 
ary condition on the surface of an inner core region. For 
non-monotonic velocity fields, the problem would need to be 
divided into monotonic flow regions and pieced together. It 
may not, however, be necessary to solve the coupled moment 
and ray equations to obtain the scattering rates. We find for 
all solutions that the linear approximation I v ~ J„ + 3-ff„/i 
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holds to high accuracy over a very broad frequency region (at 
least 100 Doppler widths for T — 10 K gas) across the line 
centre, giving the Eddington approximation value /„ ~ 1/3, 
for linear flow fields before converging to the free-streaming 
value /„ = 1 far in the wings. Allowing for a perturbed ve- 
locity flow gives close to the expected value g v = 3/5 over a 
similarly broad region. This considerably simplifies the ra- 
diative transfer computation, as solving the ray equation 
may be circumvented, requiring only solutions to the mo- 
ment equations. 

We found frequency redistribution produces solutions 
with different features across the line centre compared with 
the coherent scattering case near the source, resulting in 
variations in the scattering rate of up to ~ 50 percent about 
the coherent scattering results. Further from the source the 
solutions for a homogeneous expanding IGM were found to 
agree closely with the solution for coherent scattering, and 
roughly follow the analytically predicted r _7//3 radial depen- 
dence out to r ~ 1, beyond which the diffusion approxima- 
tion breaks down. For r > 1, the scattering rate more nearly 
approaches the free-streaming value 1 / inr 2 before becoming 
causally truncated. 

Recoils are found to suppress the scattering rate by ap- 
proximately 20 percent for a medium at a temperature of 
T = 10 K, in good agreement with estimates based on the 
diffusion approximation solution for the radiation produced 
by a uniform and isotropic distribution of sources in a uni- 
formly expanding homogeneous medium. Our computations 
extend the result beyond the diffusion approximation, and 
show that nearly the same suppression factor applies for an 
isolated source. Very near the source, frequency redistribu- 
tion modifies the suppression factor by up to 15 percent. 

In Section [5] we examined the continuum source 
Lya scattering problem allowing for inhomogeneities in the 
surrounding scattering medium, namely an overdense shell 
and a quadratic velocity profile. We found substantial devi- 
ations in the profile of the Lya scattering rate in each case 
compared with a homogeneous medium. The overdense shell 
produces not only an enhancement of the Lya scattering rate 
within the shell, but boosts the scattering rate between the 
shell and the source as well as a consequence of backscat- 
tering. A shadowed region with a deficit in the scattering 
rate compared with the homogeneous medium solution ex- 
tends beyond the shell some distance before recovering to 
the homogeneous medium value. The scattering rate pro- 
duced by the quadratic velocity profile is enhanced over the 
rate for the corresponding linear velocity profile, except near 
the outer boundary, where the rate is lower. 

As a less-contrived example of an inhomogeneous 
medium we considered a spherically symmetric density and 
velocity perturbation in a uniformly expanding homoge- 
neous IGM that satisfies the linear continuity equation. 
We found that the resulting scattering rate increases non- 
linearly with increasing values of the perturbation ampli- 
tude for perturbations beginning to become nonlinear, and 
grows exponentially with the amplitude once nonlinear as 
the velocity profile flattens near the source. We infer that 
the Lya scattering rate will depend sensitively on the veloc- 
ity structure of the IGM. 



APPENDIX A: Lya SCATTERING SOURCE 
FUNCTION FOR THE MOMENT EQUATIONS 

A particularly useful representation of the source func- 
tion for resonance-line scattering results from assuming the 
diffusion approximation or F okker-Planck approximation 
jRvbicki fc deU'Antonicl Fl994). in which a Taylor expansion 
of J v (r) under the integral in equation (133 gives 



&(r) = J„(r) + 



(A/, D ) 2 d 

2<p(y) dv 
Av v 



<p(v) 



dJJr) 



dv 



<:) 



(Al) 



<p(v) dv 



M»)Mr)] 



where we have included an additional correction term pro- 
portional to the recoil parameter e = hv ct /(\/2k^Tmc 2 ) = 
0.025(T/K) _1,/2 , arising from the effec t of atomic recoil on 
the f requency redistribution function (|Fieldl Il959bl ; iBaskd 
1 19781 ). In this section we will describe how this form of the 
source function is represented in the discretised system of 
moment equations. 

We assign a frequency grid given by the discrete set 
{vk}, k — 1, 2, 3, NF, where v\ = v max is the bluest fre- 
quency and successive values of k denote redder frequencies: 



Vl 



> V2 > Vz > 



> 1>NF- 



(A2) 



The radius grid is given by {rd}, d = 1,2,3, ...,ND, where 
n = R and tnd = Re, increasing values of the depth index 
d denote greater depths with respect to the 'surface' of the 
system at r = R: 



ri = R > r 2 > r- A > ... > r ND = Rc- 



(A3) 



In terms of the discrete grids in v and r, we write the source 
function Sk,d = S Vk (rd) as a quadrature sum: 



iJk' 



(A4) 



where Jy,d = Jv k ,(rd)- We require the coefficients 1Z k >,k,d 
for k,k' = 1,2, ...,NF and d = 1,2, ...,ND, the values of 
which are determined by a suitable discretised description 
of equation (|A1|) . 

In a system with a uniform temperature it is prefer- 
able to work with the dimensionless frequency variable 
x — (v — v a )/Avo, the offset from line centre in Doppler 
widths, and rewrite radiation quantities as e.g. J(x, r) = 
(Avn)Jv(r). Similar relations hold for higher order mo- 
ments and the source function. In these units we define 
R(x',x) = (Avd) 2 R(v',v), while the line absorption pro- 
file is expressed as <f>(x) = (Avd)ip(v). The source function 
and its diffusion approximation expansion may then be ex- 
pressed as 



S , (x, r 
+ 



(j>(x 
1 d 



2<j>(x) dx 



R(x' , x) J(x' , r) dx' = J(x, r) 



<t>{x) dJ ^ r) +2e4>(x)J(x,r) 



(A5) 



The frequency derivative terms are represented using a cen- 
tred finite difference scheme, and the coefficients multiplying 
values of Jk.d = J(xk,rd) are compared with equation (|A4|) 
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to obtain the values of lZk'.k,d 
4>k-i + 4>k+i 



k',k — 



+ 



1 - 



4>k-i 



ki^xy 

Sk'k-2 + 



k' ,k 



2cj> k Ax 



Sk',k-1 



fk{^X 
^k + l 



5y 



k',k + 2 



2<f> k Ax 



5y 



fe',fe+li 



(A6) 



where we have noted that any dependence on radius is re- 
moved for a uniform temperature medium, so that TZy , k ,d — *• 
TLy ,k- Otherwise a radial dependence must be introduced to 
the discretised values of 6. 



APPENDIX B: MONTE CARLO METHOD 

We will assume a uniform temperature medium to allow use 
of the scaled frequency offset x defined for a unique value 
of the Doppler width A^d, although the method is easily 
adapted to a variable temperature T(r) by formulating it 
in terms of the frequency variable v. We consider a pho- 
ton packet emitted/scattered with comoving frequency a; em 
at position r cm in a direction denoted by unit vector k. In 
spherical symmetry it is sufficient to designate the radius 
of emission 7"cm and the angle to the local normal # em , the 
latter having cosine /i cm = cos# em = k • r em /r cm . The path 
of the photon packet prior to the next scattering event is 
given by r = r cm + Ak, where A is the distance along the 
path of the packet. We consider a 'projected photon packet 
path' corresponding to the range A = — oo — > oo, having im- 
pact parameter r m j n defined as the minimum distance from 
the origin: r min = r cm sin6l cm = i- cm (l-^ m ) 1/2 - A spherical 
boundary of radius r > r m i n will be intersected by the pro- 
jected photon packet path at distances along the path A± 



determined from solving r • r = r with r = r c 
A± (r) = -k ■ r cm ± \r 2 - r^ m + (k ■ r cm 



Ak: 



1/2 



— /^cm^cm [T 'mi 



2 sl/2 



(Bl) 



where obviously only positive values denote possible inter- 
sections of the actual photon packet path with the sphere. 
At the radius r, the corresponding angle to the local normal 
is given by fj, = k • r/r where r = r om + A±k: 

M±( r ) = (/icmrcm + A±)/r 

= ±[l-Wr) J ] V2 ' (B2) 

At any point r along its path, the packet will have a comov- 
ing frequency corresponding to the same lab frame frequency 



x(r) 



Xem - V(r) ■ k/6 + V(rem) ■ k/6, 

■'cm fj,{r)V(r)/b + McmV(r cm )/b 



(B3) 



where fj, describes the angle to the local normal of the sphere 
at r. For the special case of a Hubble-flow velocity field, i.e. 
V(r) = Hr, this reduces to 



x(r) 



(fir — /i em rem) H/b 
\H/b 



(B4) 



where the content of the brackets simplifies following equa- 
tion (|B2f) : the frequency shift is linearly dependent on the 



distance A travelled along the path. We use these equations 
to record the frequency, path length and the angle to the 
local normal of the photon packet as it crosses a series of 
spherical boundaries of radius {rd}, the values of which com- 
prise a radial grid across which we seek the mean intensity 
J v (r) and the scattering rate P a (r). 

The packet is scattered when the accumulated optical 
depth r reaches a sufficiently large value. This optical depth 
is chosen according to an exponential probability distribu- 
tion e~ T , i.e. r = — InR for a uniform deviate R. For a 
Lya scattering opacity \ v = nuaipih 1 ), the differential opti- 
cal depth accumulated along a length dA may be expressed 
as dr = nu<y[4 l {x) / Ai^o] dA. The total optical depth is the 
integral along the path. For a static medium of constant 
temperature and neutral hydrogen density 71h, this integral 
is trivial: A = rAi^)/[nHf however, for a non-static 
medium the frequency x varies as in equation (|B3|) and no 
simple r-dependence may be obtained as a function of A. 

A less rigid approach is to assume the density, temper- 
ature and velocity of the medium are essen tially constant 
across shells, as in the Monte Carlo code of iDiikstra et all 
(2006). Under these circumstances the comoving frequency 
of the packet changes in steps as the path crosses shell 
boundaries, and the optical depth between two boundaries 
behaves as in the uniform static case. The optical depth be- 
tween boundaries at distances Xd and A<j+i along the packet 
path is 



At = : Ad 4 



A 



■Xd\ 



(B5) 



where the subscript d indicates evaluation at the intersection 
with the boundary of radius r<j, using equations (|B1|) and 
(|B3[l . The packet is permitted to progress along its path until 
the sum of the values of At exceeds the randomly selected 
total optical depth, at which stage the packet is scattered. If 
the optical depth corresponds to a total distance A ma x along 
the path of the packet, the packet is moved to a new radius 
r given by applying the cosine rule to the triangle of sides 



r and A n 



+ A 

max cm^craVax J 



\l/2 



(B6) 



Upon scattering, the photon packet is re-emitted with a new 
direction /i em and a new frequency x cm . 

We generally assume isotropic scattering, for which the 
emitted direction is obtained USlIlg f-l'cra 

= 2R — 1 for a uni- 
form random deviate R. We determined the emitted fre- 
quency from th e well-known fr equency redistribution func- 
tion Ru(x',x) (Mihalas Il978l ) using a numerical lookup 
table method. The method uses the probability distribu- 
tion p(x\x') = Rn(x' , x) / cf)(x) for the output/emitted fre- 
quency x given the input /a bsorbed f requen cy x , and is 
described in some detail in llligginsl (|2012h . An alterna- 
tive method that is more adaptable and has been widely 
adopted by other authors is based on selecting the scat- 
tering atom thermal velocity components and calculating 
the o utput frequency directly from the resulting Doppler 
shifts ijChen fc Miralda-Escudil2004r ). We found our lookup 
table method to speed up the Monte Carlo computations 
by a factor of approximately 3 compared with the direct 
Doppler shift computation method, however the required 
pre-calculation can become cumbersome. 

We estimate the mean intensity using the method of 
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iLucvl (|l999l ). based on the time-averaged number density of 
photons in a cell of a given volume. We tailor the method 
to solve problems in spherical symmetry by adopting a cell 
of volume V corresponding to the volume of the shell be- 
tween two radial grid-spheres of radius r and r + Ar, given 
by V ~ 4ivr 2 Ar. If a Monte Carlo photon packet travels a 
distance 5X in the shell between two scattering events, the 
packet occupies the shell for a time St = SX/c and con- 
tributes to the specific number density within the shell, at 
the frequency of the packet, the increment V~ 1 5t/At where 
the quantity At is the physical time represented by the entire 
Monte Carlo run. The total number density n v (r) Av comes 
from the sum of all the values of 8X between any two general 
'events', which may be a scattering event, a boundary cross- 
ing or redshifting/blueshifting between adjacent frequency 
bins, arising from packets travelling within the shell while 
having frequency within the range [v,v + Av]. The mean 
intensity is then computed from J v = [c/(47r)]n„ as 



Mr) 



ATvV(Au)(At) 



£ 5X(v,r). (B7) 



The distances 6 X between scattering events and/or 
boundary crossings are computed by following the intersec- 
tions of the path of a photon packet with the radial grid. For 
a packet that travels a short distance A ma x before scattering, 
so that the entire path is contained within the shell at r and 
the packet frequency remains within the frequency bin at v, 
we take SX — A max . However, if the path of the packet crosses 
shell boundaries or causes redshifting/blueshifting between 
adjacent frequency bins, it is necessary to separate the to- 
tal path length A max into multiple values of 6X(u, r) at the 
appropriate radii and/or frequencies. For packets crossing 
shell boundaries the SX values are obtained from differenc- 
ing the distances along the path of the packet at which the 
shell boundary intersections defined by equation (|B1|) oc- 
cur. For a packet path defined by r cm and p om the choice of 
A+ or A_ is given by either of two sequences, distinguished 
by the sign of /i cm : (i) for a packet emitted in an 'outward' 
direction with fj, em > 0, only A+ potentially represents a 
boundary crossing as A_ < and the packet may poten- 
tially cross every boundary with radius r > r cm , and (ii) 
for a packet emitted 'inwards' with /i om < 0, the solution 
A + > for all r > r m i n and A_ > for r mln < r < r cm 
and thus the potential path of the packet will cross shell 
boundaries of decreasing radii from r cm down to r m i n with A 
given by A_, then increasing radii from r m i n upwards with 
A given by A+. The values of SX are given by the differ- 
ences in sequential values of A, although we only follow the 
sequence in each case until the packet is scattered. Some ad- 
ditional calculation is needed for non-static media in order 
to treat redshifting/blueshifting between adjacent frequency 
bins while a packet is 'in transit,' i.e. in between scattering 
events and/or boundary crossings. We discuss this in appli- 
cations to particular problems within the main text. 

As an alternative to the summed path-length estimator 
for the mean intensity J v (r), we describe an estimator for 
the specific intensity I v (r, (J,). From the definition of specific 
intensity we write the differential number of photons in the 
frequency range [v, v+du] crossing the boundary at r within 
solid angle du about direction n, normal to area dA, in 
time dt as dJV = I v (r, n) dAdu dv dt. The differential solid 
angle is simply du = \d(i\ d(j> and the area dA normal to 



n is related to the differential area of the boundary dS by 
dA = dS. The total number N(r, fi, v) of photons within 
frequency range [v, v + dv] crossing the boundary of radius r 
within angular range [fi, n + dfi], is then equal to the integral 
of diV over the azimuthal angle and the boundary surface: 
N(r,n,u) = I v (r, n)\n\(f dS) (J d(j>) dfx du dt, or 

N(r,n,v) 



I v {r,n) = 



(B8) 



8ttVV|(A m )(A!/)(A£) 

where Afi and Av are the bin sizes in angle and frequency, 
and At is again the physical time represented by the Monte 
Carlo run which is related to the total number of photons 
emitted by the source. The quantity N(r, fi, v) is determined 
simply by incrementing an appropriate counter every time a 
photon packet crosses the shell at r. This method provides 
an alternative means of determining the mean intensity J„, 
which follows from angular integration, however, as it sam- 
ples only those packets that cross radial shells, we found it to 
be a less efficient means of estimating the mean intensity at 
line centre, where the path lengths are short, compared with 
the summed path-length method. This method is, however, 
required to determine higher-order angular moments of the 
intensity. 



APPENDIX C: ANALYTIC SOLUTIONS FOR A 
HOMOGENEOUS EXPANDING MEDIUM 

CI Lya Source 

In this section we review the derivation of the solution for 
a Lya point source in a homogeneous expanding medium 
jLoeb fc R vbicki 199^)- This solution provides the basis for 
the derivation of the corresponding solution for a continuum 
source. 

We rewrite equation JT]) in a dimensionless form using 
the frequency and radius variables previously described in 
Section l4~2l 



^ dr 
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dlnV 
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In the 
Ry- 



(47rf) 2 

where a = V/[H(z)r], \ = r,Xv an d V = (r*/Il)r) v . 
idealised zero-temperature limit assumed by Loeb 
bicki, atoms have zero thermal velocity and RII redistri- 
bution becomes coherent in the frame of the observer, while 
the line profile is that of natural broadening alone. Coherent 
scattering in an expanding medium means photons emitted 
at v = v a can only get redder, and thus the wing form of the 
opacity applies while the emissivity reduces to r\ v = \vJv- 
They assumed a velocity profile for the H i medium given by 
the Hubble expansion, V(r) — Hr; note that for this linear 
scaling the expression in square brackets in equation ()C1|I 
reduces to unity and 5 = 1. The zeroth and first order an- 
gular moment equations are given by: 

dJ_ 

dv 



1 d(f 2 H) 



dr 



+ 



5{v) 



6(f) 



OK 3K - 

dr f 



J dH 
dv 



(47rf) 2 



H 



(C2) 



(C3) 



A solution may be obtained in the Eddington approximation 
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K = J/3. In the diffusion limit the intensity is written as 
the Legendre expansion I v = J v + 3H v n with H v <C J„, 
corresponding to the large scattering limit f <(/, and thus 
we neglect the derivative 8H/8D in the first order angular 
moment equation. The resulting diffusion relation between 
H and J is substituted into equation ()C2|I to obtain a single 
equation for J; if we do this and change frequency variable 
to a = D 3 /9 we obtain 



dJ 
da 



1 8 



,8J_ 

Or 



5(f) 



5(D) 



(4-Kr) 2 {3a 2 ) 1 / 3 
8(f) 



(Aitf)'- 



X 8{a) 



(C4) 



where we have multiplied by the factor dD/da = (3a 2 )~ x / 3 
in the first line and used 8(D) — 5(a)/\du/da\ in the second. 
This equation is the diffusion equation in three dimensions 
with a acting as the time variable and a point source impulse 
at a = 0. The solution is given by (see Loeb & Rybicki's 
equation 21) 



3/2 



exp 



9f 
~4P 



The corresponding flux, after using Loeb & Rybicki's solu- 
tion for the mean intensity in the first order angular moment 
equation, is 



H(f,D) 
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C2 Continuum Source 



It is useful to consider the radiative transfer moment equa- 
tions applicable in the zero temperature and diffusion lim- 
its for a source term of general frequency dependence s(D). 
The combined moment equation is given by equation (|C4[) 
with 8(D) — > s(D). We use a Green's function approach to 
solve the combined moment equation by first considering 
a monochromatic source emitting at a single frequency D B 
with s(D) = 8(D - D s ) = (3a 2 ) 1/s 8(a - a s ) where a s = £ s 3 /9. 
We label the solution as G(f , a, a s ). The Green's function 
satisfies 



dG _ 1__9_ 

9cr f 2 df 
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df 
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The solution is then the 3D diffusion solution for a point 
source impulse at a — a s : 



4tt [47r(<T- CTa )] eXP [ 4 



a > a B 



G(f,a-a s ) = { ^ Y^{a-a B )\ [ 4( CT - CTs )J 

[O; cr < a s 

(C8) 

where the factor 1 / (47r) ensures unit source normalisation. 
Given G, the solution to the problem for a source of general 
frequency dependence satisfying equation (IC4|) with 8(D) — > 
s(£) is given by 



J(f,a) 



G(r,a~a s ) da s 
(3ai)y <i 



(C9) 



Thus we may construct the solution for a source of arbitrary 
frequency dependence from this result and equation (|C8|) . 

We now return to the case of a fiat source spectrum as 
an approximation to the spectrum of a continuum source in 



the IGM, for which the dimensionless source term has fre- 
quency dependence s(D) — > 1. It will be useful to specify a 
blue cutoff frequency £ m < 0, the maximum frequency (or 
minimum value of D) at which the source emits photons ca- 
pable of interacting with the Lya line, giving a step-function 
form for the source frequency distribution: 



s(a) = @(a - am) 



1 ; a > a„ 
; a < <7 n 



(CIO) 



where <r m = D m /9. We use the Green's function and the 
source frequency distribution to obtain the solution accord- 
ing to equation (|C9|) . valid for a > <r m : 



J(f,a) 



1 

3573 



1 r 

l/3( 4w )S/2 J 



a s 2,/3 G(f , a — a s ) da s 
1 



exp 



4(c 



H(f,a) 



34/3 Qf 



2(4tt) 5 / 2 



2/3/ 

a B (a 



- 2/3 G(f,a~a s )da s 



exp 



\5/2 



4(ct - a s 



da a , 
(Cll) 



da s . 
(C12) 



The forms given by equations (|49|l and (|50p are obtained 
with the change of variable u = 9f 2 /[4(D 3 — D B )]; these latter 
forms are much easier to numerically integrate. 

For v m — > — 00, equations (|C11|I and ()C12|) have the 
asymptotic series representations in t = — 4<r/f 2 



J(f,t) 

and 
H(f,t) 



2 7 \ 1/3 T(7/6) 7/3 
3 J (47r) 5 /2 r 



^ r(n + 7/6) „,,..„ 
1 + Z^ r(7/6) b ^ 2 / 3 ' w ) f 

n — 1 \ / J 



(C13) 



2/3 



r(i3/6) 10/3 2 

(4tt) 5 /2 



^ r ( 13 / 6 ) 
for \t\ <C 1, noting the general series expansion 



(C14) 



I P (f,t) = / da s a a 2/3 (ct — a s \ 

J — OO 



■p c °7[ f (° 



4V-^ /_2 \ r / +n _l\ m 



where 6(— 2/3, n) is a binomial coemcent (1, —2/3, 5/9, 
For \t\ > 1, 



(4 



r(p-i)t 



-2/3 



1 _2p- i± 
3 t 



(C16) 

where the number of terms retained in the bracketed sum 
does not exceed p. This gives the leading behaviours 

J >'^(y) I/3 W r7/3 i (C17) 
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Figure CI. Solid lines: The ratio vd log H /d log v for various f 
ranging in decades from 10~ 5 to 10 from bottom to top. Dashed 
lines: The ratio using the asymptotic expansions for |4;> 3 /9f 2 | >C 
1 and \iu 3 /9f 2 \ > 1. 



and 



H{r,t) 



(4nr) 2 



(C18) 



The diffusion approximation breaks down when dH jdv 
is no longer negligible compared with H/v 2 . The approx- 
imations equations (|C14[I and (|C18[) may be used to es- 
timate the radius at which the solution will fail. For 
|i| < 1, DdlogH /dlogu ~ 2v, while for \t\ > 1, 
59 log H/d log v ~ — (27/4) (f jv) 2 . Equating these gives as 
an approximation for the value of the frequency at which 
\vd log H/dlog v\ peaks of £ pea k ~ — (3/2)f 2//3 , with a peak 
value of \vd log H/dlog z>| pC ak — 3f 2 ^ 3 . Accordingly, the ap- 
proximation is valid only for f <C l/3 3,/2 . Comparison with 
Fig. lCll shows the peak of |z><91og H/dlog v\ exceeds 0.25 for 
r > 0.1, at which the diffusion approximation will begin to 
break down. 

The solution obtained in Section l4~4l shows that while at 
large f, \vd log H/dlog i>\ approaches unity at large values 
of \v\, for \v\ <C 1 the ratio is still small, so that dH/dv 
may continue to be neglected. Equations (|C2|I and (|C3|) then 
admit the similarity solution in t for \t\ 1 

J(r,t) = jj^m, (C19) 



(47rr) 



H{f,t) = 



3(4tt) 



j(t) + t 



*'(*) 



dt 



(C20) 



H + H 2 
2 l -t- 4 i 



if* + 



and K = J/3, where j(t) ~ 
This is found to agree well with the numerical solution ob- 
tained from the ray/moment method. 



APPENDIX D: 
SOLUTIONS 



TABLES OF TEST SUITE 



Tables of the angle-averaged intensity of the radiation field 
computed using the ray and moment method are available 



on-line for six test problems, as summarised here. Unless 
stated otherwise, the maximum frequency adopted for the 
continuum sources is x m = 1000. 

Dl Test 1: Emission line source in a static 
medium 

An analy tic solution in t he diff usion approximation is pro- 
vided by iDiikstra et al.l (|2006l ). We provide two solutions 
for RII scattering, (a) in the Eddington approximation and 
(b) solving the full ray and moment equations. The solu- 
tion using the full ray and moment equations is tabulated in 
Table [DTI The results in the Eddington approximation are 
very similar. 

D2 Test 2: Emission line source in a uniformly 
expanding homogeneous medium 

An anal ytic solution in the diff usion approximation is pro- 
vided bv lLoeb fc Rvbickil (|l999T) as well as Monte Carlo code 
results. (Also see Appendix[C]) We provide a solution for co- 
herent scattering solving the full ray and moment equations. 
The results are tabulated in Table ID2I (Vanishingly small 
values of log 10 J are indicated by a value of —100.) 

D3 Test 3: Continuum line source in a uniformly 
expanding homogeneous medium 

An analytic solution in the diffusion approximation is pro- 
vided in Appendix [C] We provide three solutions, solving 
the full ray and moment equations for (a) coherent scatter- 
ing, (b) RII redistribution, and (c) RII redistribution with 
recoils, all for a source upper frequency cutoff of x m = 1000. 
An additional solution for RII redistribution with recoil is 
also provided with x m — 1.4 x 10 s , corresponding to Ly/3. A 
temperature of T = 10 K is assumed for all solutions. The 
results are tabulated in Tables ID3ID6I 

D4 Test 4: Continuum line source in a uniformly 
expanding medium with an overdense shell 

The density profile is described by equation (|62|l . We provide 
a solution for coherent scattering solving the full ray and 
moment equations. The results are tabulated in Table [D7l 

D5 Test 5: Continuum line source in an 

expanding homogeneous medium with a 
quadratic velocity profile 

The velocity profile is given by equation (|63p . We provide 
a solution for coherent scattering solving the full ray and 
moment equations. The results are tabulated in Table |P8l 

D6 Test 6: Continuum line source in a medium 
with a self-consistent linear density and 
velocity perturbation around the source 

The density and velocity perturbations are described by 
equations (|64|l and (|65p . We provide a solution for coher- 
ent scattering solving the full ray and moment equations, 
illustrated for both a linear perturbation with Aq = 0.5 
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Table Dl. Test 1: J(r, x) (in units of 10 -9 ) for Lya source in homogeneous, static medium (T = 10 K, kq = 100, = 5): RII redistribution 



x r = 600 r = 650 r = 700 r = 750 r = 800 r = 850 r = 900 r = 950 r = 1000 



0.000 10.08642 7.96955 6.26884 4.87326 3.70389 2.70188 1.81926 

0.100 10.08641 7.96955 6.26883 4.87325 3.70389 2.70188 1.81926 

0.200 10.08642 7.96955 6.26884 4.87326 3.70389 2.70188 1.81926 

0.300 10.08641 7.96955 6.26883 4.87325 3.70389 2.70188 1.81926 

0.400 10.08642 7.96955 6.26884 4.87326 3.70389 2.70188 1.81926 

0.500 10.08641 7.96955 6.26883 4.87325 3.70389 2.70188 1.81926 



Table D2. Test 2: log 10 J(r, i>) for Lya source in uniformly expanding homogeneous medium: coherent scattering 



l°Sl0 f 


lo Sio 5 = - 1 - 5 


iogio = - 1 - 


log 10 v = -0.5 


log 10 v = 0.0 


log 10 v = 0.5 


logio = 1 -° 


iogio 5 = L5 


-3.000 


5.4131 


3.2324 


1.0482 


-1.0900 


-3.1419 


-5.1395 


-7.1773 


-2.987 


5.4119 


3.2327 


1.0483 


-1.0899 


-3.1419 


-5.1395 


-7.1773 


-2.975 


5.4103 


3.2328 


1.0484 


-1.0899 


-3.1419 


-5.1395 


-7.1773 


-2.963 


5.4085 


3.2328 


1.0484 


-1.0899 


-3.1418 


-5.1395 


-7.1773 


-2.950 


5.4064 


3.2329 


1.0485 


-1.0898 


-3.1418 


-5.1395 


-7.1773 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 



Table D3. Test 3a: log 10 J(f , x) for continuum source in uniformly expanding homogeneous medium: coherent scattering 



x logio f = -4.2 log 1() f = -3.9 log 10 r = -3.6 log 10 r = -3.3 log 10 r = -3.0 log 10 f = -2.7 log 10 f = -2.4 



-0.624 


7.53670 


6.84844 


6.15478 


5.45860 


4.76086 


4.05933 


3.35097 


-0.374 


7.53669 


6.84844 


6.15478 


5.45860 


4.76086 


4.05933 


3.35097 


-0.125 


7.53668 


6.84844 


6.15478 


5.45860 


4.76086 


4.05933 


3.35097 


0.125 


7.53668 


6.84844 


6.15478 


5.45860 


4.76086 


4.05933 


3.35097 


0.374 


7.53667 


6.84844 


6.15478 


5.45860 


4.76086 


4.05933 


3.35097 


0.624 


7.53667 


6.84844 


6.15478 


5.45860 


4.76086 


4.05933 


3.35097 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 



Table D4. Test 3b: log 10 J(f, x) for continuum source in uniformly expanding homogeneous medium: RII redistribution 



x logio? = -4.2 logio f = -3.9 log 10 f = -3.6 log 10 r = -3.3 log 10 r = -3.0 log 10 f = -2.7 log 10 f = -2.4 



-1.248 


7.58208 


6.99331 


6.36254 


5.63254 


4.82461 


4.07154 


3.36266 


-0.749 


7.58208 


6.99330 


6.36254 


5.63255 


4.82461 


4.07154 


3.36266 


-0.250 


7.58208 


6.99331 


6.36253 


5.63254 


4.82461 


4.07154 


3.36266 


0.250 


7.58208 


6.99330 


6.36254 


5.63255 


4.82461 


4.07154 


3.36266 


0.749 


7.58208 


6.99331 


6.36253 


5.63254 


4.82461 


4.07154 


3.36266 


1.248 


7.58208 


6.99329 


6.36254 


5.63254 


4.82461 


4.07154 


3.36266 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 



Table D5. Test 3c: log 10 J(f,x) for continuum source in uniformly expanding homogeneous medium: RII redistribution with recoils 



x logio r = -4.2 log 10 f = -3.9 log 10 r = -3.6 log 10 r = -3.3 log 10 f = -3.0 log 10 f = -2.7 log 10 f = -2.4 



-1.248 


7.56362 


6.94959 


6.28331 


5.51825 


4.70934 


3.97269 


3.26682 


-0.749 


7.56013 


6.94608 


6.27983 


5.51476 


4.70586 


3.96920 


3.26333 


-0.250 


7.55665 


6.94262 


6.27634 


5.51128 


4.70237 


3.96572 


3.25985 


0.250 


7.55316 


6.93911 


6.27286 


5.50780 


4.69889 


3.96223 


3.25636 


0.749 


7.54968 


6.93565 


6.26937 


5.50431 


4.69541 


3.95875 


3.25288 


1.248 


7.54619 


6.93215 


6.26589 


5.50083 


4.69192 


3.95526 


3.24939 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 
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1.00640 
1.00640 
1.00640 
1.00640 
1.00640 
1.00640 



0.00030 
0.00030 
0.00031 
0.00033 
0.00035 
0.00038 



Note: The full table is published in the el 
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Table D6. Test 3d: 
x m = 1.4 x 10 5 



i J(f, x) for continuum source in uniformly expanding homogeneous medium: RII redistribution with recoil for 



X 


log 10 f = -2.9 


logio' 1 = -2.2 log 10 f=-1.5 loj 


?10 f = -0.8 lo{ 


; 10 r = -0.1 


logio f = °- 6 logio f = L3 


-2.488 


4.46782 


2.82518 1.21451 


-0.40622 


-1.94001 


-3.42438 


-4.87538 


-1.493 


4.46080 


2.81814 1.20748 


-0.41326 


-1.94704 


-3.43141 


-4.88241 


-0.498 


4.45394 


2.81130 1.20064 


-0.42010 


-1.95388 


-3.43825 


-4.88925 


0.498 


4.44692 


2.80426 1.19359 


-0.42714 


-1.96093 


-3.44530 


-4.89630 


1.493 


4.44006 


2.79743 1.18676 


-0.43398 


-1.96776 


-3.45213 


-4.90313 


2.488 


4.43304 


2.79038 1.17971 


-0.44102 


-1.97481 


-3.45917 


-4.91018 


Note: The full table is published in the electronic version of the 


paper. A portion 


is shown here only for guidance 


regarding its form 


Table D7. Test 4: log 1() J( 


r, x) for continuum source in uniformly expanding medium with overdense shell: coherent scattering 


X 


logio f = ~ 4 - 2 


logio ' = ~ 3 - 9 !ogl0 f = - 3 - 6 lo l 


5l0 f = -3.3 lof 


r 10 f= -3.0 


logio' 1 = -2-7 


!°gio ? = -2-4 


-0.855 


7.55695 


6.94865 6.53139 


5.85581 


5.07235 


4.25632 


3.03729 


-0.513 


7.55693 


6.94865 6.53139 


5.85581 


5.07235 


4.25632 


3.03729 


-0.171 


7.55692 


6.94864 6.53139 


5.85581 


5.07235 


4.25632 


3.03729 


0.171 


7.55692 


6.94864 6.53139 


5.85581 


5.07235 


4.25632 


3.03729 


0.513 


7.55691 


6.94864 6.53139 


5.85581 


5.07235 


4.25632 


3.03729 


0.855 


7.55691 


6.94864 6.53139 


5.85581 


5.07235 


4.25632 


3.03729 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 



and extended into the non-linear regime with Ao = 2.9 (the 
latter case is not a fully self-consistent cosmological pertur- 
bation). The results are tabulated in Tables |P9l and iDTOl 
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Table D8. Test 5: log 10 J(f,x) for continuum source in homogeneous medium with quadratic velocity profile: coherent scattering 



X 


logio f = - 4 - 1 


log 10 f = -3.8 


logio r = -3.5 


logio r = -3.2 


log 10 f = -2.9 


logio f = -2-6 


log 10 f = -2.3 


-0.997 


7.80416 


7.26099 


6.64381 


5.91760 


5.09385 


4.21712 


3.30071 


-0.598 


7.80408 


7.26095 


6.64379 


5.91760 


5.09385 


4.21712 


3.30071 


-0.199 


7.80405 


7.26094 


6.64379 


5.91760 


5.09385 


4.21712 


3.30071 


0.199 


7.80403 


7.26093 


6.64378 


5.91760 


5.09385 


4.21712 


3.30071 


0.598 


7.80401 


7.26092 


6.64378 


5.91760 


5.09385 


4.21712 


3.30071 


0.997 


7.80399 


7.26091 


6.64378 


5.91759 


5.09385 


4.21712 


3.30071 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 



Table D9. Test 6a: log 10 J(f,x) for continuum source in perturbed expanding medium (Aq = 0.5): coherent scattering 



x log 10 f = -3.0 log 10 f = -2.5 log 10 r = -2.0 log 10 f = -1.5 log 10 f = -1.0 log 10 f = -0.5 log 10 f = 0.0 



-2.488 


4.86634 


3.70669 


2.54773 


1.39233 


0.21623 


-0.93978 


-2.12461 


-1.493 


4.86612 


3.70667 


2.54773 


1.39233 


0.21623 


-0.93978 


-2.12461 


-0.498 


4.86612 


3.70667 


2.54773 


1.39233 


0.21623 


-0.93978 


-2.12461 


0.498 


4.86612 


3.70667 


2.54773 


1.39233 


0.21623 


-0.93978 


-2.12461 


1.493 


4.86612 


3.70667 


2.54773 


1.39233 


0.21623 


-0.93978 


-2.12461 


2.488 


4.86612 


3.70667 


2.54773 


1.39233 


0.21623 


-0.93978 


-2.12461 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 
Zheng Z., Miralda-Escude J., 2002, ApJ, 578, 33 
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Table D10. Test 6b: log 10 J(f,x) for continuum source in perturbed expanding medium (Aq = 2.9): coherent scattering 



X 


logio f = -3-0 


logio f = "2.5 


log 10 f = -2.0 


logio r = - 1 - 5 


logio f = - 1 - 


logio f = -0.5 


logio f = 


-2.488 


5.92320 


4.76664 


3.60367 


2.43260 


1.21809 


-0.20383 


-1.97424 


-1.493 


5.92113 


4.76643 


3.60365 


2.43260 


1.21809 


-0.20383 


-1.97424 


-0.498 


5.92111 


4.76642 


3.60365 


2.43260 


1.21809 


-0.20383 


-1.97424 


0.498 


5.92111 


4.76642 


3.60365 


2.43260 


1.21809 


-0.20383 


-1.97424 


1.493 


5.92111 


4.76642 


3.60365 


2.43260 


1.21809 


-0.20383 


-1.97424 


2.488 


5.92111 


4.76642 


3.60365 


2.43260 


1.21809 


-0.20383 


-1.97424 



Note: The full table is published in the electronic version of the paper. A portion is shown here only for guidance regarding its form and content. 
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